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FOREWORD 


A  vast  amount  of  literature  exists  that  concerns  itself  with  spectrum  analysis. 

In  particular,  the  last  two  decades  have  seen  a  great  increase  of  published  theoretical 
work  as  well  as  applications  of  so-called  modern  spectral  estimation  techniques  based 
on  nontraditional  parametric  approaches.  Many  papers  dealing  with  this  subject 
make  it  difficult  even  for  the  initiated  to  make  the  connection  to  his  well  known  and 
routinely  applied  methods  based  on  the  Fourier  Transform  when  he  reads  about 
"modeling  a  random  process  to  obtain  a  parametric  spectrum”  of  the  observed  data. 

The  report  attempts  to  show  the  transition  from  the  classical  approach  of 
spectral  estimation  to  the  basic  parametric  methods.  It  also  shows  the  results  of 
these  new  techniques  when  applied  to  short  data  records  and  demonstrates  their 
capabilities  and  limitations. 


Approved  by: 


DR/y.  ETGOELLER 
Deputy  Department  Head 
Underwater  Systems  Department 
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CHAPTER  1 
INTRODUCTION 

Estimation  of  the  power  spectral  density  of  a  sigpial  is  a  popular  tool  in  signal 
analysis  for  detection  purposes  (by  enhancing  the  signal-to-noise  ratio  of  an  expected 
signal)  and  classification  (by  showing  the  distribution  of  signal  energy  across  the 
frequency  band).  Although  efforts  to  break  a  spectrum  into  its  constituting  parts  go 
back  at  least  as  far  as  Isaac  Newton  when  he  split  white  light  through  a  prism, 

Jean  Baptiste  Fourier  (1822)  established  that  any  periodic  process  can  be  interpreted 
as  an  infinite  sum  of  sine  and  cosine  terms  whose  frequencies  are  integral  multiples 
of  the  fundamental  frequency  of  the  process  (harmonic  spectral  analysis). 

It  was  Arthur  Schuster  (1897)  who  coined  the  term  "periodogram”  for  the  plot  of 
the  squared  magnitude  of  the  Fourier  coefficients  to  designate  the  spectrum.  For 
noisy  data,  however,  these  spectra  tended  to  be  quite  random  and  inconsistent,  and 
many  researchers  lost  interest  in  the  periodogram.  In  1930,  Norbert  Wiener  estab¬ 
lished  the  fundamental  relationship  between  the  autocorrelation  function  and  the 
spectrum  as  a  Fourier  Transform  pair  which  then  Robert  Blackman  and  John  Tukey 
used  as  their  approach  to  the  spectrum  (1958).  They  also  introduced  averaging  and 
windowing  to  eliminate  the  randomness  and  to  improve  the  sidelobe  structure  of  the 
original  spectra.  This  method  became  the  most  popular  spectral  estimation  tech¬ 
nique  until  Jim  Cooley  and  John  Tukey  developed  in  1965  the  Fast  Fourier  Trans¬ 
form  (FFT)  algorithm  which  calculates  Schuster’s  periodogrsun  digitally  in  a  very 
efficient  way.  This  so-called  FFT  is  now  the  most  widely  used  spectral  estimation 
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technique.  It  may  be  of  interest  here  that  Karl  Gauss  apparently  found  out  about  the 
principle  of  the  FFT  in  1805  before  Fourier  even  had  published  his  fundamental 
paper. 

Parallel  to  these  developments,  an  entirely  different  approach  was  followed  by 
statisticians,  geologists,  and  economists  who  tried  to  use  existing  test  data  to 
generate  modeled  data  with  similar  spectral  characteristics.  They  used  appropriate 
filters  whose  parameters  (parametric  analysis)  were  determined  by  the  method  of 
least-square  errors.  The  filter  is  driven  by  a  white  noise  source.  Arthur  Schuster 
tried  this  approach  as  far  back  as  1927  in  order  to  determine  the  periodicity  in  the 
cycle  of  maximum  sunspot  numbers.  This  parametric  approach  to  spectrum  analysis 
has  become  quite  popular  since  the  1960s  because  of  its  promise  of  greater  resolution 
in  transient  analysis  (John  Burg,  1967).  A  multitude  of  different  algorithms  was 
developed  since  then,  based  on  various  physical  interpretations  for  determining  the 
filter  parameters  as  well  as  on  different  mathematical  techniques  to  compute  these 
parameters  efficiently  (fast  programs). 

This  report  attempts  to  show  the  transition  from  the  classical  approach  of 
spectral  estimation  to  some  basic  parametric  methods  and  their  results  when  applied 
to  short  data  records  (transients). 

Besides  numerous  papers,  two  books  by  Steven  Kay^  and  Lawrence  Marple2 
have  served  the  author  to  familiarize  himself  with  the  material  and  to  acquire  some 
working  knowledge  for  applying  it.  Also,  the  programs  developed  by  those  authors 
have  been  used  to  calculate  the  various  signal  spectra.  These  programs  have  been 
combined  with  input  (signal  generation)  and  output  software  (display  and  file 
transfers)  so  that  they  can  be  used  as  a  flexible  interactive  package  (see  Reference  3). 
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CHAPTER 2 

REVIEW  OF  CLASSICAL  SPECTRUM  ANALYSIS 

SPECTRA  OF  ANALOG  CONTINUOUS-TIME  FUNCTIONS 


The  Fourier  series  of  a  deterministic  periodic  function  x(t)  of  period  Ti  =  1/fi  is: 


in2nf  t 

x(t)  =  ^  F(n)-e 


where  the  Fourier  coefficients  F(n)  are: 


1  f  ^1^  -jnaif.l 

F(n)  =  —•  x(t)  e  -dt  .  n  =  0,±l,±2, — 

T,  ii*-T,/2 

Because  of  tlie  periodic  nature  of  the  signal  x(t),  the  spectrum  is  a  line  spectrum 
which  exists  only  at  the  discrete  frequencies  nfi,  n  =  0,  ±  1,  ±2, . . . .  In  general,  the 
spectrum  is  a  complex  amplitude  spectrum.  As  a  special  case,  the  spectrum  becomes 
real  if  the  signal  is  an  even  function  in  time. 

The  Fourier  Transform  F(g))  of  a  deterministic  but  aperiodic  function  x(t)  is: 


1  f  * 

K(u>)  =  —  x{t) ' dl  ,  (j  =  2iif. 

2]  1  .  S  .  QC 


The  inverse  transform  reestablishes  the  original  time  function: 
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x(t)  =  F((i))  •  dw  . 

‘  CJ  S  .  <£ 


Here  the  complex  spectrum  Ftco)  is  an  amplitude  density  spectrum;  it  is  continuous  in 

CO. 


|F(u>f 


is  the  energy  density  spectrum  and 


|F(w)|^dtij 

S  *  GC; 


the  total  signal  energy. 

If  x(t)  is  not  known  exactly  at  all  times,  it  is  a  nondeterministic  (random) 
function,  and  its  Fourier  Transform  (FT)  cannot  be  determined  in  the  time  domain  (in 
theory).  Therefore,  the  key  to  the  spectrum  is  the  autocorrelation  function  ACF.  The 
Wiener-Khinchin  theorem  states:  "The  autocorrelation  function  of  a  random  func¬ 
tion  and  the  power  density  spectrum  of  the  random  function  form  an  FT  pair”  (see 
Reference  4).  The  autocorrelation  function  can  be  determined  from  the  statistical 
properties  of  the  random  function,  especially  its  first  and  second  probability  densities 
p(xi),  p(xi,  X2;  ti,  t2),  where  X)  and  X2  are  the  values  of  xft)  at  two  positions  ti  and  t2. 
If  the  process  x(t)  is  "wide-sense  stationary"  in  the  first  and  second  moment  (i.e.,  has 
constant  mean  and  variance),  the  probability  density  function  (PDF)  is  a  function  of 
the  time  difference  x  =  t2-t] .  Then  the  ACF  is  or  abbreviated  r(t): 
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r(i)  = 


x,X2-p(x,,X2;i)dXjdx2. 

..  <S 


The  parameter  t  is  also  called  the  "lag"  value  of  the  ACF.  The  ACF  deHned  by  the 
above  equation  is  an  ensemble  ACF  because  it  uses  statistics  obtained  by  observing 
all  possible  sample  functions  of  x(t)  at  time  separation  i.  This  correlation  function  is 
a  deterministic  (generally  aperiodic)  function  which  therefore  has  a  Fourier 
Transform: 


P(U))  = 


2n  j  1= 


r(i)'e 


P(a))  is  defined  as  the  FT  of  the  aperiodic  function  r(x),  and  is  therefore  the  "amplitude 
density  spectrum"  of  r(t).  From  the  definition  of  P(w)  follows  the  inverse  transform: 


r(i)  =  du» . 

‘  m  s  » <i> 


If  a  random  process  is  "ergodic,”  the  ensemble  averages  and  the  corresponding  time 
averages  are  equal.  For  this  class  of  signals,  the  time  autocorrelation  function  and 
the  ensemble  autocorrelation  f'mction  are  equal. 

The  definition  of  the  time  ACF  of  x(t)  is  given  by: 

1 

r<i)=lim  —  x(t)- xll  +  iWl 

2T  J  _t 
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when  X  =  0,  the  two  above  equations  give: 


r(0)  =  lim 
T--*® 


1 

2T 


x^Udt  = 


P(c>))du 


If  x(t)  represents  a  voltage  across  or  current  into  a  1 Q  load,  this  mean  square  value  is 
the  power  taken  by  the  load.  Therefore,  P(o))  (the  FT  of  the  ACF  of  x(t))  represents 
the  power  per  hertz  of  frequency,  or  the  power  density  spectrum  of  x(t).  Since  r(x)  is 
aperiodic,  its  Fourier  Transform,  i.e.,  the  power  density  spectrum  of  x(t),  is 
continuous  over  frequency. 


An  equivalent  definition  of  the  power  spectral  density  (PSD)  is: 


Piw)  =  E 


lim 
T— •• 


x(l)*e 


-ju>l 


dl 


-T 


where  E  means  the  "ensemble  average”  or  "expectation”  or  plain,  the  average  over 
the  terms  inside  the  curly  brackets.  This  definition  says  to  take  the  FT  of  a  given 
"realization”  xjd)  of  the  random  function  x(t)  over  an  interval  2T  as  if  it  were  an 
aperiodic  function,  square  it,  and  divide  it  by  the  record  length.  The  result  is  a 
continuous  function  in  (o  and  random.  For  any  similar  realization,  Xi(t)  of  x(t),  one 
will  get  a  different  FT  and  its  time  average  squared.  Doing  this  for  a  infinite  number 
of  realizations  of  x(t)  leads  to  an  ensemble  of  random  spectral  functions: 


2T  I 


f  T  2 

X  a)-e"^‘-dl 
-T  ' 


i  =  1,2,  . 
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Picking  a  particular  value  for  co,  the  corresponding  values  of  the  spectral  functions 
are  totally  random  and  will  not  converge  to  a  stable  spectral  value  for  increasing  time 
intervals  2T.  Therefore,  the  ensemble  average  is  taken  over  this  infinite  ensemble 
which  will  then  represent  a  statistically  stable  power  density  spectrum  as  if  it  had 
been  derived  via  the  Wiener-Khinchin  procedure.  The  proof  for  the  equivalence 
between  both  definitions  is  not  trivial  and  can  be  found  in  References  1  and  2.  The 
spectra  derived  by  both  methods  can  be  considered  "true  spectra"  compared  to 
"estimated  spectra”  which  generally  result  from  a  practical  application,  i.e.,  from  one 
realization  Xn(t)-  In  such  a  case,  one  disregards  or  only  approximates  the  averaging 
over  an  infinite  ensemble  and  obtains  the  power  spectral  estimate 


P  (to)  = 
n 


I 

2T 


X  (t)- 
n 


“  J4i)l 


■dl 


2 


which  is  a  "sample  spectrum"  (also  original  Schuster  periodogram).  The  hat  ( ‘ ) 
means  "estimate."  Although  the  mean  of  the  sample  spectrum  will  tend  to  converge 
to  that  of  the  true  PSD  in  the  limit  (T-*®),  the  variance  will  stay  constant  and 
proportional  to  the  mean  of  the  sample  spectrum.  Figures  1  and  2  depict  the 
superposition  of  the  sample  spectra  of  five  1-Vrms  white  noise  records  of  32  and  128 
data  samples.  They  show  that  the  size  of  the  statistical  fluctuations  of  the 
periodogram  remains  constant  independent  of  the  sample  size.  Only  in  the  1950s, 
when  statistical  smoothing  (ensemble  averaging)  was  applied,  the  variance  was 
reduced  resulting  in  a  more  stable  spectrum  leading  to  the  acceptance  of  the 
periodogram  approach.  The  periodogram  is  also  called  the  direct  approach  to  the 
spectrum  (i.e.,  directly  from  the  data),  whereas  the  Wiener-Khinchin,  Blackman- 
Tukey  method  via  the  ACF  is  called  indirect. 
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SPECTRA  OF  DISCRETE-TIME  FUNCTIONS 


Both  the  indirect  and  direct  methods  exist  also  if  the  signal  x(t)  is  sampled  at 
intervals  At  to  produce  a  time  series  of  samples  x(n)  for  the  integer  -«>  <  n  <  ®. 


tXi 

P(o))=  V  (indirtcl) 

k  =  -'* 


where  k-At  corresponds  to  the  lag  parameter  t  for  the  case  of  a  continuous-time 
function  x(t),  and  r(k  •  At)  is  the  value  of  the  ACF  of  x(n  •  At)  at  that  lag. 


P((0) 


E 


iim 


1 

2M  +  1 


M  I  2 


n  =  -  M 


(dirtcl) 


These  are  discrete-time,  continuous-frequency  spectra  that  can  be  shown  to  be 
periodic  in  o)  with  period  2n/At.  They  are  "true  spectra"  with  infinite  resolution  and 
an  absence  of  side  lobes  and  require,  of  course,  an  infinite  data  sequence. 


If  the  data  sequence  x(n)  is  available  only  during  a  limited  time  window  from 
n  =  0  to  n  =  N  1 ,  the  result  will  be  an  estimated  spectrum: 


(N  -  II 

P  ..  =  ^  P(k-Al) •  e”^*^*** (indirect) 

k=  _(N  -  II 


wher 


N-l-k 


f(k-At)  =  i  > 
N 


x(n)  •  x(n  +  k) 


n  =  U 
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is  called  the  biased  ACF  and 

j  N-l-k 

Wk'At)  =  -  ^  x(n)-x(n  +  k) 

.  N-k 

n=0 

is  called  the  unbiased  ACF  (k  =  0,1,2,...,N-1). 


The  direct  sample  spectrum  (estimate)  is: 


p(.,=  - 


N-l 

^  x(n)  • 


For  the  indirect  method,  a  practical  maximum  number  of  lags  k  is  always  chosen 
smaller  than  N  (suggested  by  various  authors  as  10  to  20  percent  of  N).  The  reason 
for  this  is  to  avoid  the  greater  variance  of  the  estimated  autocorrelation  associated 
with  higher  lags  k.  To  see  this,  consider  the  maximum  lag  N-l  possible  for  N  data 
points  for  which  the  ACF  estimate  is 


a.N-i) 


1 

-  -xWt-xtN-l) 


which  is  highly  variable  due  to  the  lack  of  averaging  (the  sum  degenerates  to  one 
term  only),  regardless  of  how  large  N  becomes.  On  the  other  hand,  the  maximum  lag 
should  be  chosen  large  enough  so  that  the  ACF  has  decayed  close  to  zero  and  the  ACF 
outside  this  lag  does  not  contribute  to  the  PSD. 
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A  practical  way  to  achieve  better  stability  for  larger  k  is  to  weigh  those  ACF 
estimates  less  by  multiplying  the  correlation  function,  obtained  from  the  N  data 
points,  with  a  tapered  window  function  w(k'At),  where 


and 


0  s  w(k‘At)  s  w(0)  =  1  for  k  s  knmx  =  M 

w(k-At)  =  0  for  k  >  M  ,  M  <  N-1 . 


With  this  condition, 


M 
D  1 

k=  -M 

Txr's  is  the  Blackman-Tukey  (BT)  spectral  estimate  (1958).  Many  window  functions 
exist,  but  all  are  characteristic  of  having  their  maximum  at  k  =  0  and  then 
symmetrically  tapering  off  to  a  minimum  (often  zero)  at  the  maximum  chosen  lag  M. 

The  direct  method  became  popular  after  the  Discrete  Fourier  Transform  (DFT) 
was  calculated  by  Cooley  and  Tukey  via  the  FFT  (1965).  It  is  based  on  the  discrete¬ 
time.  discrete-frequency,  Fourier-series  type  approach  for  a  periodic  signal.  The 
signal  x(n)  for  n  =  0  to  n  =  N-1  is  stored  in  memory  and  implicitly  assumed  to  be 
repetitive  with  the  fundamental  frequency 

I 

•  ~  FTai 


and  the  harmonics  fk  =  k-fi  so  that 
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u- 


u,  =  2n - 

k  N-At 


w(n-At)-*  2n  — —  •  (n-Al)  =  2ri 
N-At 


nk 

n" 


P(k)  = 


1 

N 


N-l 

^  x(n)- 

n  =  0 


2 


This  approach  leads  to  a  harmonic  analysis  of  a  random  signal  sequence  x(n),  and  the 
spectrum  is  called  a  periodogram  since  it  is  periodic  with  kmax  =  N.  This  can  be  seen 
by  calculating  the  Fourier 'Fransform  F(k  -i-N): 


1  I"'  -J2-' 

F(k  +  N)  =  -  > 


N 


>  x(n)-e 


nik  -f  N ) 
N 


N-l  .^2u- 

x(n).e 

^  n=0 

i.e.,F(k  +  N;'  =  F(k),andF(N)  =  F(0)  since  e  -*^"”  =  1.  This  is  also  equivalent  to 
saying  that  the  spectrum  is  repetitive  with  the  period  1/At. 


Smoothing  the  individual  random  sample  periodograms  is  done  by  "pseudo¬ 
ensemble  averaging."  The  given  data  sequence  x(0),  x(l), . x(N  l)  is  divided  in  P 

nonoverlapping  segments  of  Dsamples  each  so  thatP-D  =  N.  An  independent 
sample  periodogram  of  the  Pth  segment  is  then 


'  D 


D-l 


7  X  (iTt)  c 

^  p 


mk 
-j2n  — 
U 


m*0 
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The  P  independent  sample  spectra  are  then  coherently  added  to  produce  a  smoothed 
Bartlett  periodogram.  Bartlett  also  applied  a  triangular  window  to  each  data 
segment.  Welch  (1967)  then  went  one  step  further  and  permitted  data  segments  to 
overlap,  thereby  increasing  the  number  of  averaged  segments.  This  decreases  the 
PSD  estimate  variance  further,  although  the  segments  and  sample  spectra  are  not 
independent  any  longer.  He  also  applied  a  variety  of  better  windows  to  his  segments 
before  computing  tho  FFT.  This  Welch  periodogram  is  the  most  widely  calculated 
FFT  in  use  today. 

Inherent  in  the  concept  of  the  discrete-time  FI'  over  a  finite  data  sequence  is  the 
multiplication  of  an  infinite  data  sequence  by  a  uniform  (rectangular)  window. 5 
Figure  3  shows  the  development  of  the  DFT  step  by  step,  the  time  domain  on  the  left, 
the  frequency  domain  on  the  right  side.  Multiplying  the  data  x(t)  by  a  window  w(t)  in 
the  time  domain  is  equivalent  to  a  convolution  of  their  Fourier  Transforms  in  the 
frequency  domain.  Since  the  FT  of  a  rectangular  window  is  a  (sinnfT)/(nn')  function 
(Figure  3b),  the  two  lines  (positive  and  negative  frequencies)  of  a  narrowband  signal 
(Figure  3a)  are  convolved  into  two  such  (sinnlD/lnlT)  functions  resulting  in  the 
transform  X(0*W(0  where  the  symbol  ♦  signifies  the  convolution  operation  (see 
Figure  3c).  The  function  (sinnfT)/(nfr)  is  characterized  by  a  series  of  side  lobes,  i.e., 
energy  from  a  single  frequency  "leaks”  into  neighboring  frequencies  of  the  spectrum. 
The  width  of  the  main  lobe  where  it  crosses  the  frequency  axis  is  2/T  and,  at  the  3  dB 
points  down  from  the  peak,  the  width  is  about  1/T.  This  is  the  resolution  of  a  DFT  for 
a  uniform  window,  i.e.,  the  capability  to  resolve  two  frequencies  with  a  separation  Af 
=  1/T.  The  side  lobes  of  a  strong  signal  will  interfere  with  the  main  lobe  of  a  small 
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signal  if  the  frequency  separation  is  Af  =  1/T  or  less.  In  order  to  suppress  the  side 
lobes,  one  can  use  tapered  windows.  This,  however,  causes  the  main  lobe  to  widen, 
resulting  in  less  resolution. 

The  effect  of  sampling  in  the  time  domain  is  shown  in  Figures  3d  and  3e,  where 
s(t)  is  the  sampling  function  with  the  sample  pulses  At  apart  and  its  line  spectrum 
S(0  has  a  frequency  separation  of  1/At.  Figure  3e  shows  the  sampling  of  Figure  3c 
and  its  transform:  X(0*W(f)*S(f),  a  spectrum  function  still  continuous  in  frequency. 
The  step  to  the  discrete  line  spectrum  is  done  in  Figure  3f,  where  the  windowed  time 
function  in  Figure  3c  is  considered  as  one  period  of  a  repetitive  signal.  This  forces 
both  the  time  function  3£p(n)  and  its  DFT  X(fk)  to  be  discrete.  The  time  sample 
spacing  is  At,  the  frequency  line  spacing  is  Af  =  l/(N*At),  the  signal  is  periodic  in  the 
time  domain  with  T  =  N*At,  and  the  spectrum  is  periodic  with  F  —  1/At.  Both  time 
signal  and  spectrum  contain  N  terms  per  basic  cycle,  i.e.,  the  DFT  is  simply  a  one-to- 
one  mapping  ofN  terras  ofxp(n)  intoN  terms  ofX(fk). 

Figures  3d  and  3e  also  demonstrate  the  effect  of  the  sampling  rate.  The  smaller 
the  sampling  rate  1/At,  the  closer  the  signal  lobes  move  toward  each  other  until 
neighboring  lobes  start  to  overlap.  The  minimum  sampling  rate  is  the  Nyquist  rate 
of  two  samples  of  the  highest  frequency  component  in  the  signal  band  (low  pass 
band).  Below  that,  aliasing  will  occur,  i.e.,  the  appearance  of  out-of-band  high 
frequencies  as  in-band  low  frequencies  which  then  do  not  represent  the  true  signal 
anymore:  the  signal  recovered  by  low  pass  filtering  becomes  a  mo  e  and  more 
distorted  replica  of  the  original  time  function. 

The  main  lobe  structure  of  the  spectrum  calculated  for  the  DFT  with  a  uniform 
window  is  shown  in  Figure  4.  The  main  lobes  can  be  considered  to  be  a  bank  of 
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bandpass  filters  of  width  l.T  at  their  3  dB  points.  For  a  sine  wave  signal  with  a 
frequency  fa  in  Figure  4a,  which  coincides  with  one  of  the  computed  frequencies  f^, 
the  DFT  output  would  result  in  a  response  at  the  appropriate  harmonic  and  zero  at  all 
the  other  calculated  harmonics.  A  so-called  picket  fence  effect  becomes  evident  when 
the  signal  being  analyzed  lies  between  two  calculated  DFT  frequencies,  e.g.,  fbi  where 
signal  energy  appears  to  be  generated  at  the  two  neighboring  calculated  frequencies 
with  harmonic  numbers  6  and  7.  If  the  signal  is  swept  across  the  band,  then  the 
spectral  power  response  will  follow  the  curve  of  Figure  4b  showing  a  3  dB  ripple.  The 
depth  of  this  ripple  can  be  alleviated  through  zero-padding  which  generates  an 
artificial  record  length  N'  >  N  by  adding  zeroes  to  the  actual  data  record.  The 
frequency  spacing  between  calculated  harmonics  becomes  Af  =  l/(N'*At)  with  the 
result  that  the  calculated  bandpass  filters  move  closer,  resulting  in  a  spectral 
response  with  less  ripple  (see  Figures  5a  and  5b).  The  width  of  the  main  lobe  remains 
the  same  since  the  actual  data  window  is  that  of  the  original  data  record.  Zero¬ 
padding  allows  one  to  calculate  the  power  response  for  any  intermediate  frequency  by 
adding  the  proper  number  of  zeroes.  It  results  in  an  interpolation  of  the  original 
spectrum,  but  does  not  improve  the  actual  resolution.  The  spectral  resolution  is 
equal  to  the  reciprocal  of  the  signal  observation  time  or  signal  record  length  which 
determines  the  width  of  the  main  lobe  of  the  data  window  transform  (width  of  the 
bandpass  filters).  The  longer  the  observation  time  or  actual  record  length,  the 
narrower  this  filter  becomes  with  an  equivalent  improvement  in  resolution. 
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CHAPTERS 

APPLICATIONS  OF  CLASSICAL  SPECTRUM  METHODS 

Some  examples  of  classical  spectrum  analysis  follow.  They  are  based  on  a 
Fourier  Transform  of  256  points  in  the  time  and  frequency  domain.  The  number  of 
data  points  N  can  be  chosen  up  to  N  =  256.  For  N  <  256,  a  series  of  (256-N)  zeroes  is 
attached  to  the  data  record,  so  that  always  256  FT  points  cover  the  total  frequency 
band  independent  of  the  number  of  actual  data  points.  This  way  a  reasonably  smooth 
appearance  of  the  spectrum  display  is  assured.  The  sampling  of  the  frequency  axis  is 
done  relative  to  the  sampling  frequency  fr  =  f/fs,  so  that  -0.5  ^  fr  ^  0.5  where  fp  =  0.5 
represents  the  maximum  signal  frequency  possible  for  a  chosen  fg  without  aliasing 
(Nyquist  frequency).  The  vertical  axis  represents  P(f),  with  the  PSD  plotted  in  dB 
relative  to  IV'^/Hz  (PSD  =  10  logPlO/lV^/Hz).  This  is  the  display  format  of  MISA 
(see  Reference  3). 

The  periodogram  (regular  FT)  for  a  rectangular  data  window  of  width  2M  over 
N  =  32  data  points  of  a  noiseless  complex  sine  wave  of  amplitude  e  =  1 V  is  displayed 
in  Figure  6.  The  relative  frequency  is  chosen  as  fp  =  0.25.  The  results  can  be 
extrapolated  to  any  desired  frequency  within  ±0.5  fy.  The  data  record  of  N  =  32 
points  represents  N-fp  =  8  sine  wave  cycles;  the  number  of  data  points  per  cycle  is 
the  reciprocal  value  of  f,-.  The  shape  of  the  spectrum  is  the  result  of  convolving  the  FT 
of  the  rectangular  window  WR(t),  Wgif),  with  the  FT  of  the  nonwindowed  signal  s(t), 
i.e.,  S(f),  so  that  Sulf)  =  WR(f)'cS(f),  and  Pr(0  =  |SR(f)|^.  For  a  sine  wave  signal,  S(0  is 
a  Delta  function  6(f),  which  effectively  samples  the  window  function  Wr(0  to 
generate  Pr(0.  For  the  continuous  Fourier  Transform  (CFT),  WR(f)  is 
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sinnfT 

W„(0  =  CFT{w„(t)}  =  T— —  . 
**  **  itfT 


where  T  is  the  window  duration,  and  for  the  discrete  FT  it  is 


sinnf  (2M-t- 1) 

W„(n)  =  DFT{w„(n)}  =  At - - 

•'  ”  sinnf 


where  At  is  the  sample  interval  At  =  1/fs  and  (2M  + 1 )  data  points  lie  within  the 
window.  The  transition  from  the  DFT  to  the  CPT  can  be  made  by  letting  At  -»  0, 
fs  f  fs  so  that  sinnff  -►  nfp,  and  AVnfr  =  1/nfsfr  =  1/nf.  Considering  that 
T  =  (2M -1-1  )At,  sinnfr(2M-f  l)-»  sinnfT.  CFT  becomes  zero  the  first  time  for  IT  =  1, 
from  which  one  determines  the  width  of  the  main  lobe  as  2/T.  The  3  dB  bandwidth  is 
found  as  0.89/1  and  the  spectrum  peak  is  lOlogTjep  =15  dB.  The  width  of  the  side 
lobes  is  1/T.  For  the  DFl'  the  corresponding  values  are  obtained  by  letting 
T-»N  =  2M  +  1.  The  total  number  SLofside  lobes  reveals  the  length  of  the  data 
record:  SL-f2  =  N.  The  first  side  lobe,  important  for  the  detection  of  a  weaker  source 
at  a  slightly  different  frequency,  is  down  by  13.3  dB  from  the  main  lobe. 


The  effect  of  the  triangular  Bartlett  window  WB(t)  =  1  -  |n/M|  with  n  s  M  s  N/2 
can  be  seen  in  Figure  7.  The  spectral  shape  is  given  by: 


CFT{wp(t)}  = 


T  /  8in(nfr/2) 

2  \  nfr/2  / 


or 


I)FT(Wy(n)} 


sinii r  M  ,2 

r 

sm  II  f  / 
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so  that  the  main  lobe  width  4/T  is  twice  that  of  the  rectangular  window  and  the  3  dB 
bandwidth  is  1.28/T.  The  first  side  lobe  is  down  by  26.5  dB.  The  main  lobe  peaks  at 
10  logaNjep  with  a  =  0.235. 

Spectra  for  the  same  sine  wave  record  multiplied  by  a  Hanning  and  Hamming 
data  window  are  plotted  in  Figures  8  and  9,  respectively.  Both  windows  are  raised 
cosine  functions  of  the  form  (0.5  +  a)  +  (0.5-a)  *  cos|nn/M|  where  a  =  0  for  the 
Hanning  window  and  a  =  0.04  for  the  Hamming  window.  The  3  dB  bandwidth  is 
1.44/T  for  the  Hanning  window  and  1.3/T  for  the  Hanruning  window.  The  first  side 
lobes  are  down  by  31.5  dB  and  43  dB,  respectively.  The  side  lobe  ripple  across  the 
band  i  ^st  constant  for  Hamming  with  a  floor  of  32  to  36  dB,  whereas  Hanning 
provides  a  steep  roll-off  of  18  dB/Octave.  The  main  lobe  peaks  at  10  log  aNe^  with 
a  =  0.235  for  Hanning  and  a  =  0.275  for  Hamming. 

Blackman-Tukey  (BT)  indirect  spectra  can  be  seen  in  Figure  10,  based  on  the 
unbiased  correlation  function  (CF),  and  Figure  11,  derived  from  the  biased  CF,  for  a 
window  size  M  =  32  points  for  N  =  128  data  samples.  These  two  types  of  BT  spectra 
will  be  called  from  now  on  as  the  unbiased  and  biased  BT  spectra  (uBT,  bBT).  A 
triangular  Bartlett  window  was  chosen  over  the  CF.  As  seen  for  the  periodogram,  the 
window  size  M  determines  the  bandwidth  also  for  the  BT  spectrum.  However, 
whereas  the  window  size  for  the  periodogram  is  generally  chosen  as  2M  =  N  as  to 
include  the  entire  data  record,  the  window  size  over  the  CF  for  calculating  the  BT 
spectrum  can  vary  within  M  s  (N-1)  since  the  CF  has  twice  as  many  samples  as 
the  original  data  record.  The  largest  window  chosen  then  will  include  all  points  of 
the  CF. 
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A  comparison  of  Figure  10  with  Figure  6  shows  that  the  uBT  spectrum  from  64 
correlation  values  (M  =  32)  is  identical  with  the  periodogram  having  a  rectangular 
window  over  32  data  points  (N  =  32,  M  =  16).  This  is  true  for  the  uBT  with  any 
N  >  M  =  32.  Therefore,  the  uBT  spectrum  has  the  bandwidth  of  0.89yTM  for  a 
continuous-time  signal  or  0.&9/M  for  a  discrete-time  signal.  Choosing  the  biased  CF 
results  in  smoothing  the  BT  spectrum  as  seen  in  Figure  1 1.  The  side  lobe  fluctuations 
of  the  bBT  spectrum  are  calculated  as  10  log  N/M  as  can  be  verified  also  in  Figure  12 
for  N  =  32,  M  =  6.  Three  curves  are  overlaid  in  Figure  12.  One  is  the  uBT  spectrum 
for  N  =  32,  M  =  6;  the  second  is  the  corresponding  bBT  spectrum  showing  the 
reduced  negative  excursions  of  the  PSD;  and  the  third,  in  which  the  minimum 
possible  number  of  data  points  N  was  chosen  for  M  =  6  which  drives  10  log  N/M  to 
zero  and  the  bBT  spectrum  acts  as  an  envelope  follower  for  the  uBTor  the 
periodogram.  The  number  of  side  lobes  SL  across  the  frequency  band  is  directly 
related  to  M,  SL  =  M-2. 

In  all  following  plots  concerning  the  periodogram  and  the  BT  spectrum,  white 
Gaussian  noise  of  variance  on^  IV^  was  added  to  the  sine  wave  signal  of  power 
0,2  =  1V2  ,  with  a  signal-to-noise  ratio  (SNR)  of  10  log  Oip'/on^  -  0  dB.  Figure  13  gives 
a  comparison  between  a  periodogram  with  a  rectangular  window  and  one  with  a 
Bartlett  window.  Due  to  the  broader  lobe  structure,  the  Bartlett  window  smoothens 
the  periodogram  resulting  in  lower  and  wider  lobes.  This  is  true  also  for  the  signal’s 
main  lobe. 

The  same  characteristics  can  be  seen  also  for  the  Hanning  and  Hamming 
window  in  Figures  14  and  15  when  both  are  compared  to  the  Bartlett  window.  In 
particular,  the  Hamming  and  Hanning  windows  generate  almost  identical  sample 
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spectra  for  noisy  signals.  This  is  not  surprising  because  of  the  similar  structure  of 
these  two  cosine  windows. 

Figure  16  shows  overlays  of  five  periodograms  for  a  rectangular  window,  where 
each  spectrum  is  based  on  a  different  data  record  of  64  samples  to  indicate  the 
statistical  variation  inherent  in  sample  spectra  and  to  gain  some  insight  into  the 
aspects  of  improving  the  output  SNR  and  signal  detectability  with  spectral  analysis. 
The  mean  of  the  signal-plus-noise  at  fr  =  0.25  is  again  determined  by  10  log  a  Ne^ 

(a  =  1  for  WR(t)).  The  mean  of  the  spectrum  noise  mpsD  at  fr  ^  0.25  for  the 
rectangular  window  is  the  mean  of  the  input  noise  power  On^  =  1 V^,  i.e.,  0  dB.  The 
positive  noise  fluctuations  above  this  mean  are  best  described  by  the  standard 
deviation  opsD-  The  variance  is 

„  /  ,  8in2nNf  .  ^ , 

-oMl  -e  - . 

PSD  n\  \N'«in2rif^'  / 

The  envelope  of  the  second  term  within  the  parentheses  approaches  0  for  |fr|  >  0 
rapidly  with  increasing  N  so  that  practically  opsD  -» on”-  If  a  threshold  is  set  at  a 
level  equal  to  mpsn  +  opsD  =  2on^  to  define  the  equivalent  spectrum  noise,  then  this 
level  would  be  a  constant  3  dB  above  the  mean  of  the  spectral  power,  i.e.,  above  the 
mean  of  the  white  noise  input  power.  Defining  the  signal-to-noise  (S/N)  as  the  ratio  of 
the  mean  signal  to  equivalent  noise  of  the  spectrum  determines  the  SNR  as 

Ne^ 

SNR  =  101og{S/N)=  10  log  — 

2o^ 
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Since  the  SNR  of  the  input  is  10  log  (e/on)^,  the  SNR  improvement  or  gain  through 
spectral  analysis  using  the  periodogram  derived  for  a  rectangular  data  window  would 
be  10  log(N/2),  which  is  approximately  15  dB  for  N  =  64  in  Figure  16. 

Comparing  similar  plots  of  five  different  data  records  for  the  Bartlett,  Hanning 
and  Hamming  window  with  N  =  64  in  Figures  17, 18,  and  19,  it  can  be  seen  that  the 
noise  statistics  are  practically  identical  which  should  be  expected  considering  the 
spectra  in  Figures  14  and  15.  In  a  series  of  trials,  it  was  established  that  for  these 
three  windows  the  mean  of  the  spectrum  noise  mpsD  drops  by  6  dB  compared  to  the 
rectangular  window  with  a  standard  deviation  of  3  dB  above  this  mean.  However, 
since  the  mean  of  the  signal-plus-noise  peak  is  also  about  6  dB  down  (a  »  0.235, 
0.275),  the  SNR  and  the  gain  remain  practically  the  same  for  all  four  windows.  This 
result  is  not  surprising  either  since  any  window  acts  identically  on  the  noise  as  on  the 
signal.  The  major  benefit  of  windowing  with  tapered  functions  then  appears  to  be  the 
side  lobe  depression  for  signals  of  high  SNR. 

The  following  plots  are  related  to  Blackman-Tukey  spectra  for  noisy  sinusoidal 
signals  with  both  the  signal  power  and  the  noise  power  equal  to  IV^.  Figures  20  and 
21  indicate  that  no  significant  difference  exists  between  the  uBT  and  the  bBT 
spectrum  as  long  as  the  window  size  2M  ^  N.  In  Figure  20  the  data  record  is 
N  =  128,  and  in  Figure  21  it  is  N  =  64;  in  both  cases  M  =  32.  The  plots  of  Figure  22 
compare  a  uBT  and  bBT  sample  spectrum  for  the  limiting  data  record  size  N  =  33  for 
M  =  32  from  which  one  can  conclude  that  generally  the  bBT  is  a  smoothed 
representation  of  the  uBT  spectrum  for  values  of  M/N  such  that  0.5  <  M/N  <  1  with 
increased  smoothing  as  M/N  -►  1,  Figure  23  shows  that  the  uBT  spectrum  for 
N  =  M  -f-  1  is  identical  to  the  periodogram  for  N  =  M.  This  compares  directly  to  the 
signal-only  case  of  Figures  10  and  6. 
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A  high  N/M  ratio  makes  the  Blackman-Tukey  spectrum  (biased  or  unbiased) 
useful  as  a  substitute  for  ensemble  averaging  of  sample  periodogram  spectra.  In 
Figure  24  for  N  =  128,  M  =  32,  and  N/M  =  4,  this  smoothing  effect  is  equivalent  to 
averaging  four  periodograms  of  32  data  samples;  and  in  Figure  25  with  N/M  =  8,  the 
equivalent  ensemble  average  of  eight  32-point  periodograms  is  obtained.  From  an 
SNR  standpoint ,  of  course,  it  is  better  to  use  the  data  record  available  and  perform 
the  FT  operation  over  the  entire  data  set  rather  than  dividing  it  in  several  parts  and 
averaging  the  sample  spectra.  The  signal  gain  in  the  first  case  achieves  more  than 
the  noise  reduction  in  the  second  case.  However,  in  practice,  spectrum  analyzers 
(hardware  or  software)  are  designed  for  a  maximum  number  of  data  samples.  If  the 
data  record  available  exceeds  this  maximum  size,  then  averaging  is  done  profitably. 


The  results  in  Figures  26  and  27  give  an  indication  of  the  variability  of  uBT 
sample  spectra  where  five  different  data  records  were  used.  Figure  27  (N  =  64), 
Figure  26  (N  =  128),  and  Figure  25  (N  =  256)  represent  BT  spectra  (it  does  not 
matter  if  biased  or  unbiased)  for  the  same  window  size  M  =  32  with  monotonically 
decreasing  noise  empirically  established  as 


whereas  the  mean  of  the  signal  remains  constant  as  10  log  Me^  so  that 


M 

SNK  =*  10  log  — 


,2  1+2M/N 
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and  the  gain  over  the  input  SNR  is 


G  =  10  log 


M 


1-4-2M/N 


for  a  gain  of  14  dB  in  Figure  25, 13.3  dB  in  Figure  26,  and  12  dB  in  Figure  27. 


This  concludes  the  application  of  the  classical  approach  to  spectral  analysis. 
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CHAPTER 4 

THE  PARAMETRIC  APPROACH  TO  SPECTRAL  ESTIMATION 

Conventional  FT  spectral  analysis  is  based  on  a  Fourier  series  model  of  the 
data,  that  is,  the  signal  is  assumed  to  consist  of  a  set  of  harmonically  related 
sinusoids.  The  direct  approach  (FT  of  the  data)  leads  to  higher  resolution  than  the 
indirect  approach  (autocorrelation)  because  the  correlation  values  of  only  M  <  N  lags 
are  used  for  estimating  the  PSD.  Data  windowing,  either  purposely  in  order  to 
decrease  the  variance  of  a  sample  spectrum  or  unintentionally  because  the  record 
length  is  limited,  is  the  fundamental  factor  that  limits  the  spectral  resolution. 
Windowing  of  the  data  or  of  the  estimated  correlation  function  makes  the  implicit 
assumption  that  the  unobserved  data  or  ACF  values  outside  the  window  are  either 
repetitive  or  zero.  A  smeared  spectral  estimate  is  the  consequence. 

If  it  were  possible  to  "predict"  or  model  the  signal  outside  the  observation 
window  accurately,  then  the  resolution  could  be  improved.  This  is  what  the  various 
parametric  methods  attempt  to  do.  They  try  to  model  future  data  samples  based  on 
past  and  present  data  under  the  constraint  that  the  predicted  ACF  is  the  most  likely 
estimate  which  can  be  derived  from  the  observed  data.  Because  of  the  Wiener 
relationship  between  the  ACF  and  the  spectrum,  this  is  equivalent  to  modeling  the 
predicted  spectrum  as  closely  as  possible  to  the  observed  spectral  estimate. 
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Any  spectrum  can  be  obtained  from  a  white  noise  source  by  filtering  it  with  an 
appropriate  filter.  The  task  is  then  to  design  this  filter  so  that  the  output  is  a  good 
prediction  in  a  least-square  error  sense  of  the  data  to  be  expected  on  the  basis  of  past 
data. 
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CHAPTERS 

FILTERS  WITH  RATIONALTRANSFER  FUNCTIONS 

Any  linear  filter  with  a  transfer  function  H(2)  relating  the  z-transform  Y(z)  of 
the  output  y(n)  to  the  z-transform  X(z)  of  the  input  x(n)  can  be  constructed  from  the 
knowledge  of  its  poles  and  zeroes^.?  in  the  form  of  a  rational  function  in  polynomials 
of  z  (see  Appendix  A) 


H(z)  = 


N  •  z' 
Y(z)  j  =  o 


N  b  •  2  * 
1=0 


X(z)  P  ^  _ 

7  a  -z”'  a  +  /  a  -z"' 

I  O  ^  1 

1-0  1=1 


Dividing  numerator  and  denominator  by  a©  and  renaming  bi/ao  -►  bj  and  ai/ao  ->■  a;  one 
gets 


Y(z)  +  'S  a  •z”'-Y(z)=  ^  b  •z”'  X(z) 
1=1  1=0 


Each  term  z  '  •  X(z)  is  the  z-transform  of  the  input  time  series  x(n)  delayed  by  i  sample 
time  intervals,  i.e.,  x(n-i).  Equally,  z  *  •  y(z)  and  y(n-i)  are  the  z-transform  pair  of  the 
output.  Taking  the  inverse  z-transform  on  both  sides: 


P  M. 

y(n) -t-  ^a  ■y(n-i)=  ^b. -xCn-i) 
i=i  '  1=0 
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q  p 

y(n)  =  ^b.  •x(n-i)—  ^  •  y(n-i) 

i=0  i=l 

This  is  the  genc/al  input-output  relationship  of  a  digital  filter  where  the  present 
output  is  a  linear  combination  of  the  past  p  outputs  and  the  present  and  past  q  inputs. 
The  niter  coefficients  ai  and  bi  are  the  parameters  to  be  determined.  The  direct 
realization  of  the  above  expression  can  be  seen  in  Figure  28  in  which  the  block 
marked  z  ^  signifies  one  sample  delay. 

Two  special  cases  lead  to  important  filter  classes.  In  the  first  case,  all  ai  =  0 
(except ao  =  Din  which  case 


H(2)  = 


I 


1  =  0 


y{n)  =  ^  bj  •  x(n-  ii 
1  =  0 


Its  realization  is  shown  in  Figure  29.  This  is  the  digital  equivalent  of  the  classical 
"Convolution  Filter”: 


yd)  = 


I  h(i)' x(l- i)di 

J  1  =  0 


with  the  impulse  response  h(x)  replaced  by  the  coefficients  bi. 

Other  names  are  "Finite  Impulse  Response"  (FIR)  filters,  "Non-Recursive” 
filters,  'Tapped-Delay  Line”  filters,  'Transversal”  filters,  and  "Moving  Average” 
(MA)  filters.  If  all  bi  are  equal,  a  normal  average  over  a  (q  -f  D  sample  sliding  window 
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results,  otherwise  a  generalized  weighted  average.  Finally,  since  this  class  of  filters 
has  no  poles,  it  is  also  known  as  an  ”All-Zero”  filter.  Since  no  feedback  is  invo!  .,;d, 
these  filters  are  always  stable. 


In  the  second  case,  all  bi  =  0  (except  bo  =  1). 

H(z)  =  - ^ - 

p 

a  +  ^  a  -z”' 

o  ^  I 
1  =  1 

y(n)  =  -  ^  a^  •  y(n  -  i)  x(n) 

1  =  1 

Its  realization  is  shown  in  Figure  30.  This  is  called  an  “All-Pole”  filter  or  an 
“Autoregressive”  (AR)  filter.  It  is  this  type  of  filter  which  has  found  widespread 
applications  in  the  context  of  spectral  estimation  for  short  data  records.  It  models  an 
output  y(n)  based  on  the  present  input  x(n)  and  the  past  p  outputs.  Having  defined  an 
AR  filter  and  an  MA  filter,  it  follows  that,  in  the  general  case  of  Figure  28,  the  filter 
is  an  "Autoiegressive-Moving  Average”  (ARMA)  filter.  Since  for  the  AR-filters  the 
coefficient  Oq  is  just  a  gain  factor,  it  can  be  set  to  1,  and 

H(z)  =  - j - . 

1  +  Xa,z'’ 

i  s  1 


If  one  observes  .a  random  signal  s(n)  which  one  wants  to  substitute  or  model  (in 
a  spectral  sense)  by  an  equivalent  output  y(n)  with  a  modeling  filter,  then  the  input 
should  practically  be  a  random  process,  preferably  a  white  noise  process.  From  such  a 
white  noise  source,  any  desired  output  spectrum  Py(0  can  be  obtained  with  a  proper 
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filter  that  approaches  the  spectrum  PsCf)  of  the  actually  observed  signal  s(n).  From 
the  input-output  relationship  of  a  linear  filter: 

Y(z)  =  HU)X(2) 


one  gets  the  spectral  relationship: 


p  (0  =  |M(fir-p  (f)  - - 

y  *  1  p 

I  +  V 

1  =  1 


and  for  a  white  noise  input  x  =  u  and  power 


j)  (0  =  - 

y  .  p 


I  1  +  y  a 

I 


le) 


where  At  is  the  sampling  interval  in  seconds.  If  the  filter  parameters  ai  and  the  input 
noise  power  are  known,  then  the  output  power  density  can  be  determined  for  any 
desired  frequency  as  a  continuous-frequency  spectrum.  The  filter  itself  becomes  an 
equivalent  representation  of  the  spectrum  of  the  observed  signal.  The  input  noise 
power  acts  just  as  a  scale  factor. 
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CHAPTERS 

YULE-WALKER  EQUATIONS 

References  1,  2, 8, 9,  and  10  provide  the  mathematical  background  for  the 
following  material. 

In  power  spectral  estimation  the  autocorrelation  function  plays  a  fundamental 
role  (Wiener-Khinchin),  regardless  if  one  approaches  the  problem  with  one  of  the 
classical  or  parametric  methods.  It  follows  that  establishing  a  relation  between  the 
filter-parameters  (representing  the  spectrum)  and  the  ACF  will  provide  the  solution 
to  determine  the  AR  coefficients  ai.  The  development  of  this  relation  leads  to  the  so- 
called  Yule-Walker  normal  equations  (G.  Yule,  1927;  G.  Walker,  1931)  after  the  two 
statisticians  who  used  them  in  their  AR  filter  models  to  predict  trends  in  economic 
time  series.  They  were  also  derived  as  the  Wiener-Hopf  equations  to  solve  for  the 
optimum  prediction  filter. 

These  equations  derived  in  Appendix  B  are; 

r  (k>0)=  -  N  a  -r  (k-i) 
yy  '  yy 

1  =  1 


r  (k  =  0) 
yy 


p 

a  •  r  (  -i) 

*—  I  yy  o 

1  =  1 


29 


NAVSWCTR  90-236 


where  ryy(k)  is  the  value  of  the  ACF  for  a  lag  kAt  of  the  data  modeled  from  a  white 
noise  source  and  an  AR  filter  with  coefficients  ai  through  ap.  The  constraint  is  to 
devise  the  filter  such  that  the  ACF  ryy(k)  of  the  modeled  data  y(n)  be  equivalent  to 
the  ACF  rss(k)  of  the  actually  observed  data  s(n).  Therefore,  one  can  replace  ryy{k)  in 
the  expressions  above  by  rssCk)  subject  to  the  availability  of  observed  data. 

The  first  expression  sets  up  a  system  of  p  equations  rss(k  =  1)  to  rss(k  =  p)  with  p 
unknowns  ai  to  ap.  The  second  expression  provides  one  equation  with  rgslk  =0)  for 
the  only  unknown  Ou^.  In  order  to  solve  for  the  unknown  filter  coefficients  and  the 
driving  noise  source,  the  ACF  must  be  estimated  from  the  available  data.  Once  this 
has  been  done,  the  (p-1- 1)  equations  can  be  solved,  in  principle,  by  the  Gaussian 
Elimination  method;  this  requires  a  number  of  operations  proportional  to  p^,  a  rather 
time  consuming  operation.  Any  analysis  algorithm  to  be  applied  in  real-time 
requires  fast  numerical  algorithms,  and  so  the  Durbin-Levinson  Recursion  method 
was  developed  which  requires  only  p^  operations  to  solve  the  Y ule-Walker  system  of 
equations  (J.  Durbin,  1960;  N.  Levinson,  1947). 
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CHAPTER? 

LEVINSON-DURBIN  RECURSION 

The  proof  of  this  recursive  algorithm  is  lengthy  and  can  be  found  in  Reference  1 
(pp.  161-171).  Therefore,  only  the  final  equations  are  given  here  which  determine  the 
desired  AR  coefiicients.  The  number  p  of  the  coefficients  ai  through  Sp  is  called  the 
model  order.  It  signifies  the  number  of  feedback  terms  each  with  its  own  coefficient. 
The  significance  of  the  recursion  of  the  Levinson  method  lies  in  the  fact  that  it 
evaluates  the  AR  coefficients  starting  out  from  the  lowest  model  order  p  =  2  (after 
initialization  of  the  procedure  for  p  =  1),  then  working  itself  to  the  next  higher  model 
order  (p  =s  3)  and  continues  doing  so  up  to  *Se  desired  order.  At  this  point,  the 
algorithm  has  not  only  designed  the  filter  of  the  desired  order,  but  also  all  filters  of 
lower  order.  This  allows  one  to  chose  an  AR  filter  whose  model  order  appears  to  be 
most  effective  for  the  data  process  on  hand. 

Because  of  this  recursive  feature  of  the  algorithm  the  coefficients  ai  are  given 
two  subscripts:  p  for  the  model  order,  and  i  as  the  running  index  from  1  through  p. 
For  instance,  apj  =  aio,3  would  be  the  coefficient  as  of  an  AR  filter  of  model  order  10. 

The  Levinson  algorithm  starts  out  calculating  the  highest  coefficient  (i  =  p) 
with  one  formula  and  then  proceeds  to  calculate  the  coefficients  for  i  =  l,2,....,p-l 
with  a  second  formula.  A  third  expression  determines  the  power  of  the  driving  noise 
source,  also  recursively  from  lower  to  higher  model  orders. 
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The  equations  are: 


a 


pp 


P-1 


r  (p)  +  ^  a  ,  -  •  r  (p-f) 

SB  ^  p-l.C  BB 


<  =  1 


where  ^isa  dummy  parameter 


a.=a  ,+a  -a 

pi  p-lj  pp  p-l.p-1 


i=  1.2,3,  .....p-l 


a*  =  conjugate  oomplex  ofa  (in  case  of  complex  data) 


=  0- 


2 

p-l 


These  recursions  are  initialized  with 


a 


II 


and  oJ  =  (l-la,/)T^^(0) 


The  highest-index  coefTlcient  app  plays  a  significant  role  in  Burg’s  algorithm, 
described  later.  This  coefilcient  is  known  as "refK^ction  coefTicient”  kp. 


A  sample  calculation  is  given  in  Appendix  C. 
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CHAPTERS 

LINEAR  PREDICTIVE  FILTERING 

The  AR  filter  approach  as  described  above  models  the  actual  signal  s(n)  by 
filtering  white  noise  appropriately.  It  is  possible  to  give  this  modeling  approach  a 
different  physical  interpretation.  This  leads  to  the  concept  of  "Linear  Prediction”  or 
"Linear  Predictive  Coding”  (J.  Makhoul,  1975,  see  Reference  11).  In  linear 
prediction,  one  assumes  that  the  input  (white  noise  source  u(n)  in  autoregression)  is 
unknown  and  that,  therefore,  information  about  the  actual  signal  s{n)  can  be  gained 
only  from  the  past  p  outputs.  A  linear  combination  of  these  outputs  will  result  in  a 
signal  estimate  s(n)  which  only  approximates  the  true  signal  s(n): 

p 

s(n)  =  -  ')>  «  •  s(n-  i) . 

1  =  1 

Comparing  this  prediction  filter  with  Figure  30,  one  can  see  that  it  is  equivalent  to 
the  AR  filter  with  the  input  removed.  The  error  between  the  actual  value  s(n)  and 
the  prediction  s(n)  is 

e(n)  -  s(n)  -  8(n)  =  9(n)  +  ^  a^-s(n  — i). 

i  =  l 

The  coefficients  aj  are  determined  by  minimizing  the  total  squared  error  e  with 
respect  to  each  of  these  parameters,  i.e.,  by  setting  the  partial  derivative 
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The  derivation  is  given  in  Appendix  D.  The  result  is: 


—  ^  8(n)-a(n-k) ^  N  3(n-i)-8(n  — k)  =  0 


i  =  I  n 


which  leads  by  ensemble  averaging  to 


r  (k)  =  —  a  -r  (k  — i)  1  S  k 

6fi  ^  1  SS 

i-\ 


s  p 


and  the  actual  minimum  error  power  is 


r 

^  a  •  r  (i)  +  r  (0) 
nun  I  Kb  bs 


1=1 


or 


P 

r  (Oj  =  -  ^  a  •  r  (i) 

I  nun 

1=1 

These  expressions  are  formally  identical  with  the  Yule-Walker  equations  developed 
for  the  autoregressive  filter  with  the  only  difference  that  the  input  noise  power  oi,^  of 
the  AR  filter  appears  as  the  output  error  power  emin^  of  the  predictor.  The  close 
relationship  between  the  AR  model  and  the  prediction  filter  is  also  apparent  from  the 
sample  error  expression 
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p 

e(n)  =  8(n)  +  ^  a ,  •  9(n  -  i) 
i  =  l 


which  can  be  written  as 


p 

a(n)  =  -  X  *  8(n— i)  +  e(n). 
i  =  l 


This  is  the  same  equation  used  in  the  AR  process,  with  e(n)  =  u,  so  that  solving  for 
the  AR  coefficients  aj  will  determine  the  identical  parameter  set  of  the  predictor. 

The  interpretation  of  the  All-Pole  AR  filter  as  an  optimum  predictive  filter  is 
not  just  a  reformulation  of  the  same  problem  but  leads  actually  to  new  algorithms 
and  results  as  seen  later  on  for  Burg’s  method  and  others. 
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CHAPTERS 

RESOLUTION  ASPECTS  FOR  SHORT  DATA  RECORDS 

The  motivation  behind  autoregressive  spectral  estimation  or  linear  prediction 
filtering  is  to  eliminate  the  window  effect  associated  with  a  limited  data  record  that 
limits  the  resolution  of  the  analysis  and  creates  distortions  through  side  lobes.  AR 
spectral  estimation  attempts  to  extract  "enough”  i nformation  from  the  limited  data 
record  so  that  the  autocorrelation  function  can  be  estimated  beyond  the  lags 
practically  possible  for  the  given  data  record.  The  extension  of  the  ACF  is  recursive 
through  the  first  Yule-Walker  equation 


p 

r  (k)  =  -  V  a  T  (k-i) 

8H  ^  1  80 

a  s| 

I 

which  allows  one  to  calculate  the  coefficients  aj,  given  the  values  rgsfk)  for  0  s  k  s  p. 
Assuming  short-term  stationarity  of  the  signal  beyond  the  available  data  record 
these  coefficients  can  be  considered  constant,  and  extended  values  of  rgsCk)  can  be 
extrapolated  for  k  >  p.  From  an  information  standpoint,  this  modeling  approach 
permits  a  more  realistic  continuation  of  the  ACF  for  higher  lags  rather  than  setting  it 
to  zero  or  making  it  repetitive.  In  fact,  Burg  (1961)  derived  the  same  algorithm 
through  a  third  interpretation,  based  on  the  information  concept  of  "maximum 
entropy”  (see  References  10  and  12).  He  extends  the  ACF  beyond  lag  p  by  adding  the 
least  amount  of  new  information,  therefore  "maximizing  the  entropy”  of  the  process. 
This  extension  is  the  basis  for  the  higher  resolution  of  the  AR-PSD  estimation 
compared  to  the  classical  FT-based  analysis.  AR  modeling  then  is  a  process  of 
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spectrally  matching  the  observed  data  to  modeled  data,  whereby  it  should  be 
understood  that  the  modeled  time  series  data  are  generally  random  and  are  therefore 
in  no  way  matched  to  the  observed  data.  An  infinite  number  of  random  noise  time 
series  data  can  produce  the  same  ACF  and  power  spectrum. 
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CHAPTER  10 

ESTIMATION  OF  THE  AR  COEFFICIENTS 

THE  AUTOCORRELATION  AND  COVARIANCE  METHODS 

The  Yule-Walker  equations  establish  the  relationship  between  the  AR  filter 
coefUcients  and  the  autocorrelation  function  of  the  data.  Solving  them  depends  on 
knowing  the  exact  autocorrelation  values  rgsCk)  for  a  number  of  lags  k  equal  to  the 
model  order  p  of  the  AR  filter.  Generally,  this  ACF  is  not  given,  and  it  would  require 
data  over  an  infinite  time  interval.  Usually  only  a  limited  amount  of  data  is 
available  from  which  one  can  calculate  estimated  correlation  values  and  then  proceed 
to  obtain  the  coefficients  ai  using  the  Levinson  recursion.  This  is  called  the  Y ule- 
Walker  method.  Other  methods  have  been  developed  to  derive  the  filter  parameters 
from  the  data  using  the  least  square  error  criteria  of  the  prediction  filter.  As  will  be 
shown,  the  Y ule-Walker  method  can  be  thought  of  as  a  short-time  correlation  which 
depends  nut  only  on  the  lag  but  also  on  the  summation  range  n. 

If  N  data  samples  s(0),8(l),....,8(N-l)  of  a  time  series  are  known,  then  (p-»- 1)  data 
samples  are  used  to  calculate  one  error  sample  according  to 

e(ri)  =  /  a  -sCn-i)  with  a  =1. 
i=0 

Appendix  E  shows  the  system  of  linear  equations  involving  all  errors  e(n)  possible  for 
the  observed  set  of  data.  It  can  be  seen  that  the  available  data  can  be  used  to  form 
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a  maximum  of  (N  -J-  p)  error  samples:  e(0)  to  e(N  +  p-1).  However,  for  e(0),  only  one 
valid  data  sample  s(U)  is  available,  all  others  are  set  to  zero,  since  they  are  unknown. 
For  e(l)  only  the  first  two  data  samples  s(l)  and  s(0)  are  non-zero,  and  so  on.  The 
same  is  true  at  the  end  of  the  error  sample  series  where  the  last  error  sample  (N  -l-p-1) 
has  been  calculated  with  only  the  last  valid  signal  sample  s(N-l),  the  one  before  that, 
e(N-l-p-2)  with  only  two:  s(N-l),  s(N-2),  etc.  Therefore,  the  number  of  valid  signal 
samples  in  the  calculated  error  sample  series  can  be  pictured  as  in  Figure  31. 

Now  going  back  to  the  expression  which  minimizes  the  total  squared  error 


|e(n)|'* 


of  a  prediction  filter  over  a  given  range  n  of  data  samples  s(n),  one  has: 


i> 

^  ^  s(n-i)  '3(0-14)  =  -  s(n)'S(n-k)  1  s  k  s  p 

I  s  1  I)  n 


or  in  matrix  iiotation  as  found  throughout  the  literature: 


r(l.l)  r(l,2) . rd.p) 

ai 

r(l,0) 

r(2,l)  r(2,2) . r(2,p) 

32 

—  _ 

r(2,0) 

r(p,l)  r(p,2) . '■(p.p) 

ap 

r(p.O) 

where 


r(i,k)  =  ^  s(n-i)'8(n-k) 

n 
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so  that 


p 

^  'rti.k)  =  -rtk) 
isl 


If  n  goes  to  infinity,  the  sum  of  the  lagged  products  represents  the  ACF  and  the  result 
is  the  set  of  Y ule-Walker  equations.  For  a  finite  n,  however,  with  0  s  n  s  (N-1),  there 
are  four  distinct  summation  ranges  over  which  the  total  squared  error  can  be 
minimized,  as  can  be  seen  from  Figure  31.  These  four  ranges  are: 

X  or  ^  or  ^  O'"  S 
n  =  0  nsO  n-p  n  =  p 

Each  of  these  four  cases  will  lead  to  a  different  set  of  p  equations  and  result  in  a 
different  set  of  the  ai  parameters,  i.e.,  a  different  AR  filter  and  its  associated 
spectrum.  These  differences  will  become  less  and  less  significant  the  larger  the  data 
record  length  becomes  compared  to  p,  and  all  four  cases  approach  the  Yule-Walker 
method  of  the  true  ACF. 

The  first  case 


is  called  the  "windowed”  method,  since  all  sfn)  outside  the  window  0  s  n  s  N-1  have 
been  set  to  zero.  For  the  same  reason,  the  second  case 
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is  called  the  "prewindowed,"  the  third  case 


n  =  p 


the  "postwindowed,”  and  the  fourth 


n  =  p 


the  "nonwindowed”  method.  In  this  case,  none  of  the  data  for  calculating  the  error 
samples  has  been  zeroed  out.  In  much  of  the  literature  the  first  method  is  a.co  called 
the  "autocorrelation”  and  the  last  one  the  "covariance”  method.  (This  terminology 
has  old  roots,  and  semantics  is  the  reason  that  makes  reading  the  literature  harder 
than  necessary).  Of  the  four  methods,  these  two  are  found  most  frequently,  where  the 
covariance  method  leads  to  better  resolution  spectra  than  the  autocorrelation 
method.  It  is  intuitive  that  the  calculation  of  the  AR  coefficients  based  on  the  largest 
error  series  e(0)  to  e(N  -l-p-1)  is  not  optimum  because  of  the  sloped  areas  of  Figure  31 
between  0  to  p  and  (N-1)  to  (N-1  +  p)  that  are  based  on  nonexisting  data. 


THE  FORWARD-BACKWARD  PREDICTION  METHOD 

The  prediction  filter  looked  at  so  far  is  also  known  as  a  forward  predictor.  It 
estimates  future  data  from  past  data,  going  forward  in  time.  A  backward  predictor  is 
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then  a  filter  which  calculates  past  data  samples  s(n-p)  from  the  "future"  p  data 
samples  s(n-p  -I- 1)  to  s(n),  going  back  in  time.  Both  predictors  use  the  same  data 
samples,  see  Figure  32.  The  forward  prediction  is 

p 

a(n)  =  -  ^  a[-s(n-i) 
i  =  i 


with  the  error 


p 

e^n)  =  N  a|-s(n  — >). 
i-0 


Similarly  for  the  backward  prediction: 


p^ 

s(n-p)  =  -  ^  a^-s(n-p  +  i) 
1  =  1 


with  the  error 


p 

0^11)=  ^  aJ^'sCn-p  +  i), 
1=0 


For  a  stationary  process  the  forward  and  backward  AR  coefficients  aj^and  aj^  are  the 
same  (or  the  conjugate  complex  of  each  other,  aj^  =  (aj*^)*,  if  complex  data  are 
handled),  because  it  involves  the  same  statistical  information  going  forward  or 
backward  in  time.  Therefore,  one  can  combine  the  forward  and  backward  errors  in 
order  to  get  twice  the  number  of  error  samples  compared  to  the  forward-only 
non  windowed  (covariance)  method  of  range  p  s  n  s  N-1.  This  leads  to  the  so-called 
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"Modified  Covariance"  method  which  is  based  on  the  minimization  of  the  average 
combined  forward/backward  squared  error: 


p  = 


N-l  N-1 


n  =  p 


n 


P 


Minimizing  this  expression  by  differentiating  it  with  respect  to  the  aj  and  setting  it  to 
zero  yields: 


p  N-l  N-l 

—  ^  s(n-i)-s*(n  — k)  +  8(n)-s*(n-k) 

I  =  1  n  =  p  n  =  p 


+ 


p  N-l  N-l 

^  s*(n- p  +  i) *8(11  — p  +  k)  +  ^  s*(n-p)-s(n-p  +  k) 
i  =  I  n  =  p  n  =  p 


=  0 


The  expression  in  the  first  bracket  is  the  contribution  from  the  forward  error  which  is 
identical  to  that  of  the  covariance  method,  and  the  expression  in  the  second  bracket 
comes  from  the  backward  error.  This  can  be  written  similarly  as  for  the  covariance 
method: 


p 

^  •  r(i,k)  =  -r(k)  I  s  k  s  p 

1  =  1 


where 


N-l 

r(i,k)  -  N  s(n  — i)*8*(ii  — k)  +  3*(n  — p  +  i)'  sin  — p  +  k)j 
n  =  p 


and 
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N -1 

r<k)  =  ^  ^s{n)*s*(n  — k) -h  8*(n-p)-s(n-p  +  k)^ 

II  =  p 


THE  BURG  METHOD— MINIMIZING  THE  REFLECTION  COEFFICIENTS 


This  method  was  derived  by  John  Burg  (1967)  and  has  become  a  widely  used 
technique  to  determine  the  AR  coefficients.  It  should  not  be  mistaken  for  Burg's 
maximum-entropy  interpretation  of  AR  spectral  estimation,  mentioned  earlier.  It  is 
strongly  related  to  the  forward-backward  prediction  by  minimizing  the  combined 
forward/backward  squared  errors.  However,  while  those  errors  ep^(n)  and  epHn)  were 
calculated  straight  forward  for  a  given  model  order  p,  in  Burg’s  method  they  are 
recursively  derived  from  the  errors  ep./(n)  and  ep-i^n)  of  model  order  (p-1)  by 
requiring  that  the  coefficients  Spj  follow  the  Levinson  recursion 


=  a  ,  -I-  k  • 

p-i,i  p 


p-  I.P- 


wherek|  =  app  is  the  "reflection  coefficient”  as  mentioned  earlier.  In  Appendix  E, 
the  expressions  for  the  recursive  errors  are  derived: 

c*^(n)  =  (n)  +  k  ,(n-)) 

p  p-l  p  p-i 

e^n)  =  e*”  (n-l)  +  k  (n) 

p  p-l  p  p-l 


The  average  combined  squared  error  for  the  nonwindowed  case  is  then 
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E 

p 


N-1 


,(n)+k  -e**  (n-1) 

p-1  p  p-1 


2 


+ 


e**  (n-l)  +  k  (n) 

p-1  p  P-1 


which  is  only  a  function  of  kp,  since  the  lower  order  prediction  error  powers  ep-/(n) 
and  ep.i^(n)  are 'known  (having  been  calculated  starting  from  the  first  order  model). 
Therefore, 


aE  aE 

p  ^  p 

a(Re(k  ))  '*a(Jin(k  )) 

p  p 


allows  one  to  solve  for  kp  as  the  only  unknown.  This  results  in: 


N  -  1 


-2  y  e*"  ,(n)  (e‘’  .(n-l))* 

—  p-i  \  p-1  / 


k  = 


n  =  p 


P  N-1 


;  2 


N-1 


n=p  n=p 


(Proof,  References  1 ,2) 


The  initialization  of  the  Burg  method  is  done  by 

) 

=  r  (0)  =  —  |x(n)l'‘^ 

O  68  N1  ’ 

n  =  rt 

e*^(n)  =  B(n)  n  =  1 ,2,3,..  ,(N  —  1 ) 

O 

e‘’(n)  =  8(n)  n  =  0,1 .2,..  ..,(N  -  2) 

O 


Based  on  these  initialization  errors  the  value  of  ki  can  be  calculated.  In  order  to 

determine  k2,  all  errors  ei^(n)  and  ei^(n-l)  for  n  =  p,p  +  l . ,N-1  must  be  calculated 

from  the  recursive  error  expressions 
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e^n)  -  ,(n) -h  k  ,(n-l)  n  =  p-l- 1  1 

p  p->  p  p-i 

and 

eSn)  =  (n-l)-Pk  ,(n)  n  =  p,p-H,....,N-2 
p  p->  p  p-i 

With  k2  determined,  the  errors  can  be  undated  again  to  give  the  error  series  for  e2^(n) 
and  e2^(n-l)  which  are  used  to  compute  ka,  and  soon,  until  kp  of  the  desired  model 
order  p  has  been  reached.  After  each  of  these  cycles,  the  Levinson  recursion  is  used  to 
find  the  coefficients  api  from  the  k’s: 


and  the  minimized  error  power  output  is 


a 


2 

P 


Finally,  the  spectral  estimate  is 


P  (0  = 

c 


O 


2 

P 


1 


P 


i^l 


2 


]  -I- 
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This  completes  the  Burg  algorithm.  It  should  be  mentioned  briefly  here  that 
the  two  error  recursion  formulae  lead  to  a  diir?rent  interpretation  of  the  prediction 
error  filter,  the  so-called  lattice  filter,  depic^d  in  Figure  33.  The  input  signal  s(n?  is 
the  observed  signal;  the  error  signals  ep^(n)  and  ep^(u)  are  generated  as  output.  The 
lattice  filter  is  optimally  matched  to  the  input  signal  in  the  sense  that  the  output 
error  tends  to  go  to  zero  if  the  filter  is  perfectly  matched  to  the  input  spectrum,  i.e., 
the  spectral  estimate  as  calculated  above,  is  equal  to  the  input  spectrum.  For  a  signal 
with  slowly  changing  statistics,  the  coefficients  a;  can  be  made  to  change,  in  effect 
tracking  the  input  spectrum,  such  that  the  output  is  always  minimized.  The  lattice 
filter  has  found  its  use  in  adaptive  filtering. 
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CHAPTER  11 

APPLICATIONS  OF  PARAMETRIC  SPECTRUM  METHODS 

RESOLUTION  CAPABILITY  IN  THE  (NEARLY)  NOISELESS  CASE 

The  following  is  rela  ,ed  to  narrowband  signals  (ideally  sinusoidal)  in  white 
noise.  As  stated  earlier,  the  record  length  T  of  an  observed  signal  inherently  puts  a 
limit  on  the  spectral  resolution,  i.e.,  on  the  capability  to  resolve  two  closely  spaced 
frequencies.  This  limit  of  resolution  is  generally  taken  as 


For  a  given  signal  frequency  f  and  a  record  length  T  the  number  C  of  cycles  in 

Tis 


C  =  Tf 


A  higher  frequency  (f  +  Af)  will  allow  a  larger  number  of  cycles  (C  +  AC'  in  the 
same  record  T 


c  -f-  nc  =  T  (f+  no  i.e.,  nr  =  nc/r 

This  incremental  frequency  Af  is  equivalent  to  the  limit  of  resolution  for 
AC  =  1  cycle.  This  means  two  frequencies  can  be  resolved  if  the  higher  frequency 
signal  contains  at  least  one  more  cycle  than  the  lower  frequency  signal  in  the 
record  T. 
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Figure  34  and  35  depict  the  real  part  and  the  magnitude  of  a  complex  transient 
sinusoid  containing  two  components  of  relative  frequency  0.25  and  0.3125  of  equal 
amplitude  and  with  a  relative  phase  shift  of  160  degrees  for  a  record  length  of 
N  =  16  data  points.  This  results  in  four  cycles  at  the  lower  and  five  cycles  at  the 
higher  frequency  and  the  one-cycle  difference  necessary  for  spectrally  resolving  the 
two  components. 

The  periodogram  is  shown  in  Figure  36,  and  i  t  indicates  the  importance  of  the 
relative  phase  between  components  for  short  duration  signals.  For  a  phase  shift  of 
335  degrees,  the  two  frequencies  are  clearly  resolved,  but  biased  compared  to  the  true 
frequencies  indicated  by  the  two  dashed  vertical  lines.  At  270  degrees,  both 
components  start  to  separate,  still  biased,  and  at  160  degrees  just  one  spectral  peak  is 
located  at  the  average  frequency.  The  periodogram  method  clearly  requires  more 
than  a  one-cycle  difference  to  reliably  resolve  two  frequencies  of  difference  Af,  or 
stated  differently,  the  record  length  must  be  larger  than  1/Af. 

The  Burg  method  has  been  chosen  for  calculating  the  parametric  spectrum  of 
Figure  37  for  the  most  difficult  phase  shift  of  160  degrees,  showing  the  superiority  of 
the  autoregressive  modelling  approach. 

In  the  following  three  illustrations,  Figures  38, 39,  and  40,  the  resolution 
capability  of  three  spectral  techniques,  the  FT  based,  the  Burg,  and  the  Modified 
Covariance  method,  is  determined  more  accurately.  In  each  case,  the  lower  relative 
frequency  is  0.25.  The  upper  one  is  varied  in  order  to  accommodate  a  desired  cycle 
difference  AC.  The  two  amplitudes  are  equal,  and  some  noise  is  added  to  prevent 
singularities  (SNR  =  26  dB). 
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Figure  38  is  a  plotof  the  periodogram  for  three  cases,  where  AC  is  1.3, 1.4,  and 
1.5  cycles,  corresponding  to  a  frequency  of  0.33125, 0.3375,  and  0.34375.  The  most 
unfavorable  relative  phase  with  respect  to  resolution  was  chosen.  The  PSD  is  plotted 
on  a  linear  rather  than  on  the  ordinary  dB-scale  which  allows  better  recognition  of 
the  two  PSD  peaks  to  develop.  It  can  be  seen  that  a  minimum  difference  of  1.4  to 
1.5  cycles  is  necessary  for  an  FT  to  resolve  two  frequencies. 

Figure  39  shows  the  Burg  spectrum  of  order  p  =  4  for  AC  =  0.5  and 
AC  =  0.7  cycles,  corresponding  to  f  =  0.28125  and  f  =  0.29375.  An  excess  of 
0.7  cycles  will  separate  the  two  frequencies,  however,  with  a  bias  depending  on  the 
relative  phase  A4>.  It  can  be  shown  that  this  bias  varies  as  a  function  of  sin  A<P. 

Figure  40  is  a  plot  of  the  Modified  Covariance  spectrum  for  AC  =  0.3  and 
AC  =  0.4  cycles,  corresponding  to  f  =  0.26875  and  f  =  0.275.  This  method  then 
requires  an  excess  of  only  0.4  cycles  to  resolve  two  frequencies.  Also,  a  significantly 
smaller  frequency  bias  than  for  the  Burg  method  exists  relative  to  phase  over  the 
entire  2n-range. 


RESOLUTION  OF  SMALL  NOISY  SIGNALS  IN  THE  PRESENCE  OF  LARGE 
SIGNALS 

Under  noisy  conditions,  the  ideal  resolution  capability  demonstrated  above  will 
be  degraded.  Given  are  N  =  16  data  points  of  a  sinusoid  of  relative  frequency 
f]  =0.25(Ni  =  4  cycles)  and  1-volt  amplitude.  A  second  sinusoid  of  the  same  sample 
size  is  chosen  such  that  its  frequency  f2  coincides  with  the  frequency  of  the  first 
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sidelobe  maximum  of  the  periodogram  at  f  =  0.34  (N2  =  5.44  cycles).  Its  amplitude  is 
0.1  volt.  White  noise  is  added  with  a  standard  deviation  of  0.1  Vrms,  so  that  SNR(fi) 

=  20  dB  and  SNR(f2)  =  0  dB. 

Figure  41  is  a  representation  of  the  FT  based  power  spectrum  for  five 
independent  data  records  indicating  the  effect  of  the  additive  noise.  No  clue  is  found 
in  these  spectra  about  the  presence  of  the  second  sinusoid. 

Figure  42  shows  the  same  five  simulations  using  a  Hamming  window.  Again  it 
will  be  quite  diHlcult  to  identify  and  locate  the  second  small  signal  component.  In 
Figure  43  the  Burg  algorithm  with  a  model  order  p  =  8  was  applied  to  the  same  five 
data  records  with  the  result  that  in  four  trials  the  small  signal  can  be  identified  using 
simple  thresholding  and  judging  the  consistency  of  peaking  in  the  spectrum. 

A  similar  experiment  was  done  selecting  the  frequency  f2  =  0.405 
(N2  =  6.5  cycles)  to  coincide  with  the  second  sidelobe  maximum  of  the  large  signal  at 
fi  =  0.25.  The  periodogram  with  a  rectangular  window  was  calculated  and  plotted  in 
Figure  44;  the  Hamming  window  was  applied  in  Figure  45.  Not  much  information 
can  be  gained  from  both  sets  of  spectra  about  the  presence  of  a  second  signal.  The 
Burg  algorithm  of  model  order  p  =  8,  however,  provides  a  spectrum  in  Figure  46  in 
which  both  frequencies  can  be  identified  and  located. 

In  a  third  experiment,  the  second  signal  was  hidden  in  the  minimum  between 
the  main  lobe  and  the  first  side  lobe  of  the  FFT-based  spectrum  of  Figure  47  at 
f2  =  0.3125  (N2  =  5  cycles).  Figure  48  shows  the  spectrum  using  again  the  Hamming 
window.  The  Burg  spectrum  (p  =  8)  is  displayed  in  Figure  49  with  a  similar  result  as 
for  the  two  previous  trials. 
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These  experiments  allow  the  conclusion  that  small  noisy  signals  are  difTicult  to 
detect  across  an  extended  part  of  the  spectrum  and  particularly  so  close  to  the 
mainlobe  of  a  large  signal  using  the  conventional  Fourier  method.  This  gap  is  filled 
by  the  parametric  method  of  AR  spectral  analysis. 


COMPARISON  OF  PARAMETRIC  METHODS 

The  previous  series  of  tests  was  expanded  to  include  the  other  parametric 
algorithms  discussed  in  this  report,  the  autocorrelation,  the  covariance,  and  the 
Modified  Covariance  method,  together  with  the  Burg  algorithm.  The  first  relative 
frequency  is  again  at  0.25, 16  data  points  form  a  record  of  four  cycles;  for  the  second 
signal  5.5  cycles  were  selected  leading  to  f2  =  0.344,  well  within  the  resolution  limits 
of  parametric  spectra  for  any  phase  shift  between  the  two  components.  The  SNR  of 
the  first  signal  is  +  20  dB,  that  of  the  second  is  0  dB. 

Figure  50  shows  what  the  autocorrelation  algorithm  can  achieve  at  best.  Even 
at  these  relatively  high  model  numbers  p  =  8  and  10,  this  method  proves 
unsatisfactory  to  represent  a  realistic  spectrum. 

The  results  of  the  covariance  method  are  plotted  in  Figure  51  for  p  =  4, 5, 6.  At 
p  =  4  the  resolution  is  too  small,  and  the  spectrum  peaks  only  at  the  main  signal 
frequency.  For  the  next  higher  order  model,  the  second  component  starts  to  show  up, 
and  it  is  well  established  for  p  =  6.  In  either  case,  the  second  peak  occurs  somewhat 
off  the  actual  frequency.  Also  note  that  for  p  =  6  the  weak  component  is  calculated  to 
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have  a  higher  spectral  value  than  the  main  component,  indicating  that  the  model 
order  is  too  lar^e.  Going  to  still  higher  orders  will  cause  the  spectrum  to  break  up 
into  more  and  more  peaks  unrelated  to  the  actual  signals. 

The  Modified  Covariance  method  was  used  for  Figure  52.  Also  for  this 
algorithm,  model  order  4  generates  only  the  main  peak;  p  =  5  indicates  the  second 
sinusoid,  and  p  =  6  produces  a  satisfactory  spectrum.  The  Burg  algorithm  used  in 
Figure  53  produces  similar  spectra  for  the  same  model  orders.  These  two  algorithms 
have  been  found  to  be  equally  effective  for  various  signal  and  noise  conditions;  both 
are  relatively  robust  in  the  sense  that  they  often  (not  always)  tolerate  a  higher  model 
order  without  breaking  up.  The  Modified  Covariance  method  has  an  edge  over  Burg 
with  a  somewhat  better  resolution  capability,  and  it  is  also  less  affected  by  phase 
variations.  However,  the  Burg  algorithm  is  faster.  For  a  data  record  N  =  64  it  took 
5  seconds  on  a  PC  to  calculate  the  spectrum  for  p  =  24,  but  20  seconds  for  the 
equivalent  Modified  Covariance  model. 


SPECTRAL  PERFORMANCE  AT  LOW  SIGNAL-TO-NOISE  RATIOS 

Generally  the  benefit  of  high  resolution  is  derived  from  AR  spectral  analysis  at 
high  SNR  conditions.  The  following  five  sets  of  spectra  in  Figures  54  through  58  are 
the  result  of  five  simulations  of  a  16-sampIe  signal  at  f  =  0.25  and  of  white  noise 
with  an  SNR  =  0  dB.  In  each  set,  the  model  order  is  p  =  6.  The  autocorrelation 
(Figure  55),  the  Modified  Covariance  (Figure  57),  and  the  Burg  method  (Figure  58) 
show  a  detection  performance  comparable  to  the  periodogram  (Figure  54),  however. 
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with  better  resolution.  The  covariance  method  of  Figure  56  yields  highly  unstable 
spectra  that  are  very  sensitive  to  the  individual  noise  record. 

Going  to  still  smaller  signals  with  an  SNR  =  -6  dB  for  the  same  N  =  16  data 
record,  no  useful  information  is  obtained  any  longer  from  any  spectral  method, 
parametric  or  Fourier  based.  Figures  59  through  63  show  the  results.  Coherent 
averaging  of  the  five  spectra  of  each  set  will  improve  the  performance;  however,  this 
would  imply  longer  data  records  which  may  or  may  not  exist.  This  concludes  the 
testing  of  the  AR  spectral  estimation  methods  discussed  in  the  report  and  their 
comparative  performance. 
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CHAPTER  12 
CONCLUSIONS 

Whenever  new  methods  are  tried  and  added  to  the  pool  of  existing  signal 
processing  tools,  it  is  desirable  to  be  aware  of  their  capabilities  and  limitations.  What 
do  they  allow  us  to  do  that  could  not  be  accomplished  before  with  the  tools  available? 
This  is  true  also  when  considering  the  so-called  "modern  spectral  analysis"  compared 
to  "conventional  Fourier  analysis.” 

FFT-based  spectral  analysis  is  the  most  robust  method  for  widely  ranging 
signal  and  noise  conditions.  Particularly  for  low  to  moderate  SNR ,,  the  periodogram 
with  its  numerous  window  options  can  be  well  adapted  to  the  given  signal  and  noise 
characteristics  and  will  give  satisfactory  results.  Also,  it  does  not  make  any  demands 
on  the  bandwidth  of  the  signal  and  performs  well  for  broadband  or  narrowband 
signals.  The  analysis  of  a  broadband  process  is  done  best  through  the  conventional 
FFT  approach.  The  identification  of  narrowband  components  is  possible  even  for 
transients  as  long  as  the  cycle  excess  between  comp' rents  is  at  least  a  few  cycles, 
theoretically  down  to  one  cycle.  All  of  this  can  be  done  in  real  time  because  of  the  fast 
FFT  algorithm. 

When  all  is  said,  a  few  well-dellned  problems  remain  which  cannot  be  solved 
conventionally  because  of  properties  inherent  to  the  FT  when  applied  to  time  limited 
data.  The  Fourier  Transform  is  strictly  defined  for  unlimited  data  and  is  therefore  a 
mathematical  rather  than  an  engineering  concept  When  dealt,  with  r.  a  mathe¬ 
matical  context,  there  are  no  resolution  limits  or  side  lobe  effects.  However,  all 
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engineering  applications  involve  data  that  exist  only  during  a  given  observation 
time,  or  that  are  intentionally  truncated,  i.e.,  the  data  are  windowed.  A  time-limited 
window  function  is  transformed  by  the  FT  into  an  unlimited  frequency  space  through 
more  or  less  pronounced  side  lobes.  If  somewhere  in  these  side  lobes  another  signal  is 
present  smaller  than  the  side  lobe  level,  it  cannot  be  detected. 

The  modeling  approach  of  the  autoregressive  analysis  is  capable  of  solving  such 
a  problem.  By  substituting  the  observed  data  by  a  model  of  the  process  that 
generated  those  data,  one  has  access  to  an  unlimited  source  of  equivalent  data.  This 
model  is  derived  from  the  information  of  the  actual  signal,  essentially  its  associated 
correlation  function.  Since  no  limits  are  imposed  on  the  output  of  the  model,  the 
spectral  estimation  can  be  done  without  generating  side  lobes.  A  peak  in  the 
spectrum  related  to  a  sinusoid  in  the  signal  will  decay  asymptotically  to  zero,  so  that 
another  smaller  sinusoid,  even  in  the  vicinity  of  the  main  peak,  can  still  be  detected. 

The  model  is  actually  a  filter  which  is  designed  such  that,  for  each  sinusoid 
contained  in  the  signal,  a  pole  at  the  corresponding  frequency  is  generated.  This 
gives  rise  to  a  resonance  effect  that  is  very  sensitive  to  frequency  and  that  generates 
a  strongly  pointed  peak  in  the  transfer  function  resulting  in  better  resolution.  This 
high  resolution  makes  it  possible  to  separate  two  components  in  frequency  even  if 
their  cycle  difference  is  less  than  one  cycle,  the  theoretical  limit  for  conventional 
analysis. 

These  are  the  two  essential  characteristics  of  parametric  AR  spectral  analysis 
that  provide  a  new  tool  in  signal  processing  and  fill  a  true  gap  existing  in  classical 
analysis.  They  also  indicate  the  problems  associated  with  this  approach.  First,  it  is 
strictly  line  oriented  due  to  poles  in  the  transfer  function.  Second,  one  has  to  make 
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assumptions  about  the  number  of  lines  in  the  s  sectrum  corresponding  to  the  model 
order,  i.e.,  one  has  to  have  apriori  information  Uiat  is  not  normally  available. 

Several  criteria  exist  for  choosing  the  proper  model  order.  No  attempt  has  been 
made  in  this  report  to  deal  with  model  order  selection  since  this  is  a  subject  by  itself. 

It  may  suffice  here  to  state  that  the  existing  criteria  are  not  hard  and  fast;  they 
rather  provide  some  guidelines.  Generally,  parametric  spectra  do  not  appear  to  be 
very  sensitive  to  model  order  in  large  SNR  cases  (the  main  working  domain  for  AR 
analysis)  if  one  is  interested  primarily  in  frequency  identification;  resolution, 
however,  improves  with  increasing  order.  Needless  to  say  that  computation  time 
balloons  for  higher  model  order  interfering  with  many  real  time  applications. 

Finally,  it  should  be  pointed  out  that  a  peak  in  the  AR  spectrum  resulting  from 
a  sine  wave  component  is  not  proportional  to  the  power  of  the  sinusoid  as  it  is  for  the 
periodogram,  rather  it  is  proportional  to  the  square  of  the  power.  The  integrated  area 
under  the  peak,  however,  still  represents  correctly  the  power  as  it  should  be  for  a 
power  spectral  density  curve. 
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FIGURK  1.  FIVK  PEKIODOGRAMS  OF  1  VKMS  WHITE  NOISE.  32  DATA  SAMPLES 
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DEVELOPMENTOF THE  DISCRETE  FOURIER  TRANSFORM  X  (fKiOF 
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(A)  DISCRfcTE  FOURIER  TRANSFORM  AS  A  SET  OF  BANDPASS  FILTERS 
I  POWER  RESPONSE 
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(B)  RESULTING  PICKET  FENCE  EFFECT 


FIGURE  4  DISCRETE  FOURIER  TRANSFORM  RESPONSE 


fk 

(A)  CLOSER  SPACING  OF  FILTERS  THROUGH  ZERO  PADDING 
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(B)  REDUCTION  OF  PICKET  FENCE  EFFECT 


FIGURE  5  DISCRETE  FOURIEKTRANSFORM  RESPONSE  WITH  ZERO  PADDING 
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FIGURE  6  PERIODOGRAM  OF  A  COMPLEX  SINUSOID,  32  SAMPLES.  UNIFORM  WINDOW 
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FIGURE  10.  BI,ACKMAN  TU  KEY  SPECTRUM,  UNBIASED,  128  DATA  SAMPLES,  64 
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FIGURE  II.  BLACKMAN-TUKEV  SPECTRUM.  BIASED.  128  DATA  SAMPLES. 64  SAMPLE  BARTLETT  W  INDOW 
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nCURE  12.  UNBIASED  AND  BIASED  Bi.ACKMAN-TUKEY  SPECTRA.  12  SAMPLE  BARTLETT  WINDOW 
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FIGURE  13.  PERIODOGRAMOFASINUSOID  AND  WHITE  NOISE.  SNR  =  OdB.  32 DATA  SAMPLES, 
BARTLETT  VS  UNIFORM  WINDOW 
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FIGURE  14.  PERIODOGRAM  OF  A  SINUSOID  AND  WHITE  NOISE.  SNR  =  0  dB.  32  DATA  SAMPLES. 
BARTLETT  VS  HANNING  WINDOW 
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BARTLETT  VS  HAMMING  WINDOW 
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FIGURE  20  BLACKMAN  TUKEY  SPECTRA  OF  A  SINUSOID  AND  WHITE  NOISE.  SNR  =  0  dB, 
128  DATA  SAMPLES,  G4  SAMPLE  BARTLETT  WINDOW 
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FIGURE  21  BLACKMAN-TUKEY  SPECTRA  OF  A  SINUSOID  AND  WHITE  NOISE,  SNR  =  OdB, 
64  DATA  SAMPLES,  64  SAMPLE  BARTLETT  WINDOW 
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FIGURE  22  BLACKMAN  TUKEY  SPECTRA  OF  A  SINUSOID  AND  WHITE  NOISE,  SNR  OdB, 
33  DATA  SAMPLES,  64  SAMPLE  BARTLETT  WINDOW 
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FIGURE  23  PERIODOGRAM  (32  DATA  SAMPLES)  WITH  UNIFORM  WINDOW  AND  BLACKMAN  TUKEY  SPECTRUM 
(33  DATA  SAMPLES)  WITH  64  SAMPLE  BARTLETT  WINDOW.  SNR  =  OdB 
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FIGURK  24  PKRIODOGRAM(32  DATA  SAMPLES)  WITH  UNIFORM  WINDOW  AND  BLACKMAN-TUKEY  SPECTRUM 
(128  DATA  SAMPLES)  WITH  64  SAMPLE  BARTLETT  WINDOW,  SNR  =  OdB 
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FIGURE  25  PERIODOGRAM  (32  DATA  SAMPLES)  WITH  UNIFORM  WINDOW  AND  BLACKMAN  TUKEY  SPECTRUM 
(256  DATA  SAMPLES)  WITH  64  SAMPLE  BARTLETT  WINDOW,  SNR  =  0  dB 


FIGURE  27  BLACKMAN  TUKEY  SPECTRA,  5  RECORDS.  64  DATA  SAMPLES  EACH,  WITH  64  SAMPLE 
BARTLETT  WINDOW,  SNR  =  OdB 
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FIGURE  34.  REAL  PART  OF  A  COMPLEX  TRANSIENT  («6  DATA  SAMPLES)  COMPOSED  OFTHE  FREQUENCIES 
0.25  AND0.3125OF  EQUAL  AMPLITUDE  AND  leO-DEGREE  PHASE  SHIFT 
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FIGt'RE35.  MAGNITUDE  OF  A  COMPLEX  TRANSIENT(16DATASAMPLESlCOMPOSEDOFTHEFREQUENClES 
0.25AND0.3I25OFEQUAL  AMPLITUDE  AND  IfiO-DEGREE  PHASE  SHIFT 
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FIGURK  37  BURG  SPKCTRUM  (P  =  2)  OK  TRANSIENT  OF  FIGURES  34, 35 
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FIGURE  38.  PERIODOGRAMOFTRANSIENTOTH  OF  EQUAL  AMPLITUOKMB  DA  I  A  SAMPLES)  COMPOSED  OF  TWO 
FREQUENCIES  FI  =  0.25  AND  F2.  BOTH  OK  EQU  AL  AMPLITUDE 
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BOTH  OF  EQUAL  AMPLITUDE 
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FI  =  0  25  AND  F2,  BOTH  OF  EQUAL  AMPLITUDE 
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FiGURH43  BURG  SPECTRA  (P  =  8)  OF  TWO  SINUSOIDS  (FI  =  0  25.  F2  =  0  34)  AND  WHITE  NOISE,  5  RECORDS, 
16  SAMPLES  EACH,  SNR  (FI)  =  20dB,SNR(F2)  =  OdB 
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FIGUKh:44  PKRfOI)()GRAMOFTWOSINUSOII)S(Fl  =  0  25.  F2  =  0  405)  AND  WHITK  NOISE.  5  RECORDS, 
Ifi  SAMFM-KS  EACH,  UNIFORM  WINDOW.  SNR(Fl)  =  20cIR.  SNR  (F2)  =  OdB 
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FIGURK45.  PERIODOGRAM  OF  TWO  SINUSOIDS  (FI  =  0.25,  F2  =  0.405)  AND  WHITE  NOISE.  5  RECORDS, 
16  SAMPLES  EACH,  HAMMING  WINDOW,  SNR  (FI)  =  20dB,SNR(F2)  ==  OdB 
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FIGURE  46.  BURG  SPECTRA  (P-8)  OF  TWO  SINUS0IDS(F1  =  025.F2=  0  405)  AND  WHITE  NOISE,  5  RECORDS. 
16SAMPLRSEACH,SNR(F1)  =  20dB,SNR(F2)  =  OdB 
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FIGURE  48  PERIODOGRAM  OF  TWO  SINUSOIDS  (FI  =  0  25.  F2  =  0  312)  AND  WHITE  NOISE.  5  RECORDS. 
16SAMPLES  EACH.  HAMMING  WINDOW.SNR(Fl)  =  20dB.SNR(F2)  =  OdH 
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FIGURE  49  BURG  SPECTRA 
16  SAMPLES  EA 
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r'GURESO.  AUTOCORRELATION  SPECTRA  OF  TWO  SINUSOIDS  (FI  =  0.25,  F2  =  0.34)  AND  WHITE  NOISE, 
16  SAMPLES,  SNR  (FI)  =  20dB,SNR(F2)  =  OdB 
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KIGURE51  COVARIANCE  SPECTRA  OFTWO SINUSOIDS  (FI  =0  25,  F2  =  0.34)  AND  WHITE  NOISE, 
leSAMIM.ES,  SNR(Fl)  =  20dB.SNR(F2)=0dB 


2.40 
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KIGURK.‘)2  MODIKIKI)  COVARIANCE  SPECTRA  OFTWO  SINUSOIDS  (FI  =  0  25,  F2  =  0  34)  AND  WHITE  NOISE, 
16SAMPI.ES,SNR(Fi)  =  20dB,SNR(F2)=OdB 
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KIGURK53  BURGSPECTRAOFTWOSINUSOII)S(Fl  =0  25,F2  =  0  34)  ANDWMITK  NOISK, 
1 6  SAMPLES.  SNR  (Fi )  =  20  dB,  S’ ,  R  (F2)  =  0  dB 
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FIGUKK54  PKRIODOGKAMOFSINUSOUXF  =0  25)  ANDWIIITKNOISK.S  lli:COEU)S, 
1 6  SAM IM-KS  FACE! ,U N EFOE<M  WIN E)OW ,  SNEt  =  0  dH 
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FIGURE  55  AUTOCORRELATION  SPECTRA  (P  =  61  OF  A  SINUSOID  (F  =  0  25)  AND  WHITE  NOISE, 
5  RECORDS.  16  SAMPLES  EACH,  SNR  =  0  dR 
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FIGURK  56  COVARIANCK  Sl’KCTRA  =  8)  OF  A  SINUSOID  (F  =  0  25)  AND  WHITK  NOISK,  5  RKCORDS, 
1 6  SA  M  PDFS  K  ACl  I ,  SN  R  =  0  dB 
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SKKCORDS,  lCSAMPLKSKACH.SNR  =  OdH 


xIOi 
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KKIUKK  58  BURG  SPECTRA  (P  =  6)  OF  A  SINUSOID (F  =  0  25)  AND  WHITE  NOISE,  5  RECORDS,  IfiSAMl’EES  EACH,  SNR  =  OdB 
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UNIPORM  Window,  SNR  = -6  dR 
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5RKCOR[)S,  16  SAMPLES  EACH.  SNR  = -6  dB 
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FIGURE  61 .  COVARIANCE  SPECTRA  (P  =  6)  OF  A  SINUSOID  (F  =  0.25)  AND  WHITE  NOISE,  5  RECORDS. 
16  SAMPLES  EACH, SNK  =  -6dB 


IGURt:  62  MODIFIEDCOVARIANCK  SPECTRA  (P=--6)  OF  A  SINUSOID  (F  =  0.25)  AND  WHITE  NOISE, 
5  RECORDS.  16  SAMPLES  EACH,  SNR=  6  dB 
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KIGUKK  63  BURG  SPKCTRA  (P^  6)  OF  A  SINUSOID{F  =  0,25)  AND  WHITE  NOISE,  5  RECORDS.  16  SAMPLES  EACH.  SNR=  6  dR 
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APPENDIX  A 

POLE-ZERO  FORMAT  OF  FILTER  TRANSFER  FUNCTION 


This  is  to  show  the  equivalence  of  the  transfer  function  H(z)  as  a  ratio  of 
polynomials  in  z  as  used  in  the  report: 


N  z  ' 

H(z)  =  -  (A-1) 

p 


and  as  a  ratio  of  products  in  which  factors  of  the  form  (z  -  z[)  and  (z  -  pi)  represent 
zeroes  and  poles  of  the  transfer  function.  From  Equation  (A-1)  one  gets 


HU)  = 


^  +  b,z-'  bjZ-*  -h  b^z-^  -h 


a  +  a  z  *  +  a„z~^  u„z~^ 

0  1  2  3 


■i-  b  z'** 
q 


+  a 


IfO<q<p,  p-q  =  N  >0,  then 


b  z^+b-Z**  *-fb,.  z**^ 


H(z)  = 


P  P— •  A  P-2  . 

a  z*^  a,  z*^  +  a„  z'^ 


+  b 

4 


a 


b  z**  +  b,  z**  *  +  b„z'*~^  +  •  •  ■  -t-  I) 


a  z^  4  a.z**'*  4  a„z’’~^  4  •  •  •  4  a 
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b 

N  0 


+  -3 


+  _I  jq-i  +  J  2^-2  ^ 
b  b 

u  o 


’P  -H  —  zP“'  —  zP"^  + 


z**  +  B,  z**-’  -e  B-z''-^  + - +  B  ,  z  -*-  B 

,  Q  ,  _ ! _ 2 _ q-J _ q 

zP  A  zP~*  +  A,zP“^+  •  •  ■  -h  A  z  -h  A 

1  2  1  p 


where  the  numerator  and  denominator  are  ordinary  polynomials  in  z.  This  can  be 


written  as 


H(z)  =  G  •  z' 


N  — 


(Z-P,)  (Z-Pj)  (Z-Pg) 


(z-z  ) 
q 

(z-p  ) 
p 


=  G  •  z 


n 

N  ,  _ 

p 

n 

i  =  > 


where  G  is  a  gain  factor  and  implies  N  zeroes  at  the  origin  of  the  z-  plane. 


Ifp<q,  q-p  =  M>0,  then 


fl 


'  11 


with  M  poles  at  the  origin  of  the  z-  plane. 


A-2 
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APPENmXB 

DERIVATION  OF  THE  YULE-WALKER  EQUATIONS 

The  autoregressive  (AR)  filter  models  the  present  output  y(n)  through  a  linear 
combination  of  the  past  p  outputs  y(n-i)  for  i  =  1,  2, 3,  p,  and  the  present  input 

x(n)  where  x(n)  is  a  white  noise  source. 


p 

y(n)  =  —  ^  •  y(n-i)  f  x(n) 

i  =  I 


Multiplying  by  y(n-k)  gives 


y(n)  •  y(n-k)  =  —  ^  a.  •  y(n-i)  •  y  (n-k)  +  x(n)  •  y(n-k) 

I  ::  1 


Ensemble  averaging  E{  }  on  both  sides  and  interchanging  the  E-operator  and 
E-operator  (ergodic  process  assumed)  results  in: 


NAVSWCTR  90-236 


where  ryy  (k)  is  the  output  autocorrelation  function,  and  rxy  (k)  is  the  input-output 
cross  correlation.  From  the  fundamental  input-output  relationship  in  the  time 
domain  by  which  the  input  is  convoluted  with  the  impulse  response 

QD 

y(n-k)  =  ^  h(i)  ■  x(n-k-i) 

is  — » 


one  gets 


(k)  =  E  |x(n)  •  y(n-k) 


=  E 


N  x(n)  •  h(i)  ■  x(n  — k  — i) 


=  y  h(i)-E 


I  a: 


)  (B 

c(n)  •  x(n-k  — i)  =  ^ 


h  (i)  •  r  (k  +  i) 
xx 


r  (k  +  i)  =  •  6(k  +  i)  , 

XX  t 

since  the  autocorrelation  function  (ACF)  of  white  noise  x(n)  is  the  delta  function  of 
strength  Therefore, 


r  (k)  =  ^  h(i)  •  6(k  -I-  i) 

ty  1  .i — 

I  s  -a> 


8  (k  -h  i)  as  a  function  of  i  is  equal  to  zero  for  all  i  except  i  =  -k,  for  which 

-  k  +  t 

'y  8(k  i)  =  1  ,  withi;-»0  , 

i  =  -  k  -  c 
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so  that 

-k  +  c 

r  (k)  =  •  h(-k)-  y  6(k  +  i)  =  0^  ■  h(-k) 

jty  X  -to  X 

i  =  -k  -c 

The  value  for  hl-k)  can  be  found  from  the  basic  definition  of  the  z-transform  H(z) 
for  h(n): 

H(z)  =  y  h(rk)  •  2~"  =  y  h(n)  •  z~" 
n  =  -  “  n  =  0 

since  fora  realizable, causal  filter  h(n<0)  =  0  ,  i.e.,  h(-k)  =  0  for  k  >  0  . 

H(z)  =  h(0)  •  1  +  h(l)  .  z"'  +  h(2)  •  -k . 

lim  H(z)  -  h(0) 

Also: 


l  +  a,  •  z  '  -k  •  z  ^  • 

SO  that 

h(0)  -  1,  r  (k)  =  0^  for  k  =  0 

ry  I 

and 

h(-k)  =  0,  r  (k)  =  0  fork  >  0 

*y 


liin  H(z)  = 


lim 

I  -•  ® 


p 

1  +  y 

i  =  1 


=  lim 
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and  finally 


r 

r  (k)  =  —  ^  a  •  r  (k— i)  k  =  1 ,2,  •  •  •  •  p 

>7  >  yy  .  .  H 


r  (k)  = 
yy 


—  ^  a  •  r  (k-i)  +  for  k  =  0 

^  1  yy  x 


This  is  the  set  of  Y ule-Walker  equations,  the  first  of  which  results  in  a  group  of  p 
eauations  to  solve  for  the  p  unknowns  ai,  a2,  ,  ap  given  the  values  of  the  desired 

ACF.  The  second  expression  for  k  =  0  is  just  one  equation  to  solve  for 

Example  of  equation  for  k  =  1; 


ryyd)  =  -  +  a2-r^(-l)+  a3T^(-2)  +  •..•  +  apT^^d-p) 


which  can  also  be  written  as 


^  (0)  +  a-  r  (-l)  +  a,-r  (-2)  +  --  --  +  a-  r  (l~p)-0 

yy  lyy  2yy  3yy  pyy*^ 


and  similarly  for  k  =  2, 3,  ,  p. 
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The  entire  set  of  (p  +  1)  equations  can  be  written  in  the  short  form  of  matrices 
and  vectors  as 


This  is  the  form  in  which  the  Yule-Walker  equations  are  frequently  found  in  the 
literature.  One  can  see  that  the  terms  of  the  diagonals  in  the  ryy-matrix  are 
identical.  This  indicates  a  so-called  Toeplitz  structure  of  the  matrix  which  allows  one 
to  solve  such  a  system  of  equations  with  a  method  known  as  the  Levinson-Durbin 
recursive  algorithm. 


B-5/B-6 
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APPENDIX  C 

EXAMPLE  OF  THE  LEVINSON  RECURSION  FOR  P  =  4 


In  order  to  calculate  the  four  coefHcients  a4i ,  a42,  a43,  and  a44  of  the  fourth  order 
AR  filter,  one  starts  out  (after  initialization  which  gives  an  and  o]^)  with  order 

p  =  2: 


®22 


r  (2)  +  a,,-  r  (1) 

K&  II  SS 


a  =  a  +  a  *8 
21  11  21  11 


(>- 

\ 

22 

2 


'  2 

I  ■  °1 


These  are  all  parameters  of  model  order  2.  Now  proceeding  to  model  order  p  =  3: 

r  (3)  +  (a  •  r 

bB  \  21  6 


33 


a_  +  . 

a 

21 

33 

+ 

a 

22 

33 

a 

\ 

33 

Model  order  p=‘4  (with  known  parameters  for  p  =  3): 


C-1 
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These  last  five  parameters  completely  describe  the  desired  AR  filter  to  model  the 
data.  The  autocorrelation  values  rgg  (k)  must  be  known. 


C-2 
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APPENDIX  D 

DERIVATION  OF  THE  LEAST  SQUARE  ERROR 
OK  THE  PREDICTION  FILTER 

The  error  e(n)  between  the  actual  signal  sample  s(n)  and  its  prediction  s(n)  is 

H 

e(n)  =  s(n)  -  s(n)  =  s(n)  +  ^  a  s(n-i) 

1  =  1 

The  totaled  squared  error  over  a  given  time  interval  of  n  samples  is 

^  e^(n)  =  ^  ^s(n)  +  ^  a,  '  8(0-1)^ 
n  n  1  =  1 

=  ^  (n)  +  2  ^  s(n)  ^  a.  •  s(n-i)  +  ^  ^  a,  •  8(n-i)j 

n  n  i  =  1  n  i  =  1 

In  order  to  find  the  parameter  set  {ai}  which  minimizes  this  expression,  one  calculates 
the  partial  derivatives  with  respect  to  all  a’s  and  sets  them  to  zero.  In  order  not  to 
confuse  the  derivative  index  and  the  delay  index  i,  one  chooses  another  index  for 
calculating  the  derivative,  say  k.  The  derivative  of  the  first  term  is  zero,  since  it  is 
independent  of  all  a.  The  derivative  of  the  second  term  can  be  written  as 
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The  third  term  has  the  form  z  =  y^(x)  for  which 


dz 

dx 


dz 


dy  dy 

—  =  2y  •  — 
dx  dx 


The  derivative  of  this  term  car  then  be  written  as 


n  i  =  I  '^“k'k  =  l  ' 


where 


d 


aj-s(n-l)  +  *"''-(-a|^-s(n-k)  +  ••• 


=  s(n  — k) 


SO  that  the  third  term  is 


p 

I  a.  •  8(n-l)  •  8(n  — k) 

n  1  =  1 


for  k  =  1 , 2,  •  •  ■  • 


The  total  derivative  of  the  squared  error  is  then  set  to  zero: 


p 

2  ^  8(n)  •  8(n  — k)  +  2  V  ^  •  sin  — i)  •  8(n-k)  = 

n  n  I  «  1 


D-2 


•  8p'  8(n-p) 
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The  ensemble  average  of  this  expression  leads  to  the  requirement  for  the  minimum 
mean  square  error: 


E 


^  8(n)  •  8(n-k)j  =  -  ^  “i 

n  1=1 


E 


s(n-i)  •  8(n-k) 

n 


which  is 


p 

r  (k)  =  -  a  •  r  (k  -  i)  for  k  =  1 , 2,  3,  •  •  •  •  p 

as  —  1  bb  »  i  »  r 

I  =  1 
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APPENDIX  E 

SET  OF  ERROR  EQUATIONS  FOR  N  DATA  SAMPLES 


Expanding  the  error  equation  for  a  prediction  filter  of  order  p 


p 

e(p.)  =  s(n)  +  ^  •  s(n  — i)  =  ^  a.  s(n-i) 

i  =  1  i  =  0 


one  gets  the  set  of  (N  -I-  p)  error  samples  possible  for  an  observed  set  of  N  data,  s(0) 
through  s(N-l).  All  unknown  data  for  (N-1)  <  n  <  0  are  set  to  zero. 


e(0)  =  aos(O)  +  ai  +  a2rf^  +  . 

e(l)  =  aos(l)  +  ai  s(0)  +  a2»^*B  + . 

e(2)  =  aos(2)  -I-  ais(l)  +  a2s(0)  +  . 

1 

1 

..  -1-  ap-i  *^4rp)  apH-p) 

.  -4-  ap.i  +  apS-H-p) 

.  ap-is-tS-qj)  +  ap«-fS-i^ 

1 

1 

1 

1 

e  (p-1)  =  aos(p-l)  +  ai  s(p-2)  a2s(p-3)  . . 

1 

1 

.  +  ap.]  s(0)  -1-  ap9-(-4^ 

e{p)  =  aos(p)  -»■  aj  s(p-l)  a2S(p-2)  +  .. 

e(p-»-l)  =  aos(p-hl)  ai  s(p)  +  a2S(p-l)  +  . 

.  +  ap-i  s(l)  -1-  ap8(0) 

.  +  ap-i  8(2)  +  Bpsd) 

1  1 

I  Covariance  Method  ! 

1  1 

e  (N-2)  =  ao  s(N-2)-hai  s(N-3)-t-a2s(N-4)  -f-  . . 
e  (N-1)  =  aos(N-l)-t-ai  8(N-2)-l-a2s(N-3)  +  . . 

. +  ap-i  s(N-l-p)  +  apS(N-2-p) 

. +  ap-i  s(N-p)  +  ap8(N-l-p) 

e (N)  —  aosfW  +  ai  e(N-1) -»■  a2  8(N-2) ... 

e(N-t-l)  =  ao«fN-»-4^-J-ai  8W)  +  a2s(N-l)-l-  . 

e  (N  +  2)  =  ao9(N  I  g)-f-ai  s(N  <  l)-Ha2»fN)-f 

1 

1 

- +  Bp-i  s(N-p-l- 1)  +  apS(N-p) 

. .  ap-i  s(N-p-4-2)-t-apS(N-p  +  l) 

. . .  +  ap-i  s(N-p-»-3)  +  ap8(N-p-l-2) 

1 

:  : 

e  (N  +  p-3)  =  ao8(N  I  p’3)-t-ai  b(N  1  p-4)-h . +  ap-i  8(N-2)-t-  apS(N-3) 

e  (N -H p-2)  =  ao 8(N-<-p-2)  + . -f  ap.2  8^N)+ ap-i  8(N-1)  +  apS(N-2) 

e  (N-t-p-1)  —  ao 8(N  1- p*l)-t-  . +  ap-2 a(N  <  1)  +  ap  i  S(N)  +  ap8(N-l) 

Autocorrelation  Method 

E-l/E-2 
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APPENDIX  F 

DERIVATION  OF  BURG’S  ERROR  RECURSION 


The  errors  of  the  forward  and  backward  predictors  are  (with  app  =  kp) 


p-i 


E*^(n)  =  8(n)  +  V  a  •  8(n-i)  =  s(n)  +  %  a  .  •  s(n-i)  +  k  •  8(n-p) 

P  —  pi  ^  pi  p 


I  =  I 


=  I 


P.  * 

^(n)  =  s{n  — p)  ^  — p  +  i) 


I  =  1 


Substitution  of  Spi  by  Levinson’s  recursion 


a  =  a  ,  .  +  k  •  a  , 

pi  p-i,i  p  p-i,p-i 


gives 


ej(n)  =  8(n)  +  ^  (“p-l.i  +  ‘‘p  ‘ 

I  =•  1 


*  ) 

p-l.p-i/ 


s(n  — i)  +  k  •  8(n  — p) 
P 


p_i  ,  p-\  ^ 

=  B(n)  +  ^  a  ,  •  8(n-i)  +  k  I  a  ,  •  s(n  — U  +  8(n  — p) 

^  p-1,1  p\  ^ 

1=1  1*1 


where 


p-i 

B(n)  +  V  a  ,  .  •  8(n-i)  =  e'  ,  (n) 

P  —  l»  I  P-  I 

i  *  1 


FI 
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Now  substituting  p-i  =  k,  n-i  =  n-(p-k) 
and  with  is=l-»k  =  p-l,  i=(p-l)-»k  =  l 

1  ^ 

-►  ^  s(n-p  -f  k)  8{n-p) 

k  =  p-I 

Renaming  k  i  and  noting  that 

I  p-i 

1=1. 

k=p-I  k=l 


* 

>  a  <  a(n-i)  +  8(n-p) 

^  p-i,  p-i 

i  =  J 


this  expression  becomes 


p->  ^ 

(n-p)  4  X  ®  -1  i  *  8(n-p  +  i) 


i  =  l 


=  s^(n-l)-(p-l)j+  ^  Bpj  .  •  8^(n- 1 )  -  (p- 1 )  + 


SO  that  the  forward  error  is  recursively 

e*^  (n)  =  e*^  <n)  +  k  •  e’’  ,  (n  -  1 ) 

p  p-l  p  p-i 

The  derivation  of  the  recursion  backward  error  is  similar  and  results  in 

e^n)  =  ,(n-l)  +  k*  ,  (n) 

p  p-l  p  p-l 


F-2 
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